Estimating WEPP Cropland Erodibility Values From Soil Properties

William J. Elliot; Dennis C. Flanagan

doi:10.13031/ja.15218

Article Request Page ASABE Journal Article
Estimating WEPP Cropland Erodibility Values From Soil Properties
William J. Elliot^1,*,†, Dennis C. Flanagan²
Published in Journal of the ASABE 66(2): 329-351 (doi: 10.13031/ja.15218). 2023 American Society of Agricultural and Biological Engineers.
¹ USDA-Forest Service, Rocky Mountain Research Station, Moscow, Idaho, USA.
²USDA-Agricultural Research Service, National Soil Erosion Research Laboratory, West Lafayette, Indiana, USA..
^*Correspondence: welliot@moscow.com
^†Retired.
Submitted for review on 1 June 2022 as manuscript number NRES 15218; approved for publication as a Research Article and as part of the Soil Erosion Research Symposium Collection by Associate Editor Prof. Anita Thompson and Community Editor Dr. Kyle Mankin of the Natural Resources & Environmental Systems Community of ASABE on 15 December 2022.
Mention of company or trade names is for description only and does not imply endorsement by the USDA. The USDA is an equal opportunity provider and employer.
Highlights
Very fine sand is the single most important variable for predicting rill and interrill erodibility values.
Slope steepness is the most important property for predicting critical shear for rill erosion.
Erodibility prediction equations with the best goodness-of-fits have both soil texture and mineralogy terms.
Equations to predict erodibility from texture, organic carbon, cation exchange, slope, and taxonomy are proposed.
Abstract. In the late 1980s, the USDA Agricultural Research Service, along with other federal agencies and multiple universities, collaborated to develop a new physically based soil erosion model, the Water Erosion Prediction Project (WEPP) Model. The WEPP model was intended to replace the Universal Soil Loss Equation and was to include estimates for upland runoff and erosion, sediment delivery to first order channels, and runoff and sediment routing through a downstream channel network. The WEPP technology estimated erosion from raindrop splash and sheet flow (interrill erosion) and concentrated channel flow (rill erosion). To make these erosion estimates, WEPP required new soil erodibility values for interrill erodibility (Ki). rill erodibility (Kr), and critical shear (??c) for concentrated flow erosion. The WEPP Core team determined that they needed to estimate these three erodibility values from measurable soil properties for a wide range of soil conditions. To develop relationships between WEPP soil erodibility variables and other soil properties, a field study was carried out using rainfall and runoff simulation to measure the three erodibility values for 36 soils. Sites were identified on croplands from Washington to Georgia and Maine to California, USA for erodibility measurement. Concurrently, the USDA Soil Conservation Service (SCS) carried out detailed soil surveys and laboratory analyses for all sites to provide a large database of soil physical, chemical, and engineering properties. Correlation and regression analyses were carried out to develop relationships between SCS measurable soil properties and WEPP soil erodibility values. This article provides a summary of the field procedures, data analyses, and subsequent predictive equations that were developed. The predictive equations that were finalized in the WEPP User Summary used sand, very fine sand, clay, and organic carbon contents to predict cropland soil erodibility, but the Coefficient of Determination (r2) values were 0.55 or less. More complex predictive equations were developed with soil physical, chemical, mineralogical, and geomorphic properties, with r2 values up to 0.81. Most of the better predictive equations included terms for soil texture and clay mineralogy, often with additional chemical properties. A set of simplified erodibility equations using only the readily available properties of soil texture, organic carbon, cation exchange capacity, slope steepness, and taxonomic order were derived for use within the WEPP Model, with r2 values greater than 0.5 for all three equations for estimating soil erodibility from measurable soil properties.
Keywords. Critical Shear, Interrill Erodibility, Rill Erodibility, Soil Properties, WEPP.
A corrigendum to this article can be found at: Elliot, W. J., & Flanagan, D. C. (2024). Corrigendum to “Estimating WEPP cropland erodibility values from soil properties.” J. ASABE, 67(5), 1241-1244. https://doi.org/10.13031/ja.16069
Soil erodibility is a qualitative or quantitative property that describes the ability of a given soil to be detached and transported by an erosive agent, such as wind or water (Ellison, 1947; Lal and Elliot, 1994). Soil erosion by water is due to interactions among climate, soil properties, topography, and both natural and human disturbances (Foster and Lane, 1987; Huffman et al., 2013; Wischmeier and Smith, 1978). Natural and human activities and time have resulted in a wide range of soils with physical and chemical properties that contribute to soil erodibility (Buol et al., 2003; Huffman et al., 2013; USDA Natural Resource Conservation Service, 2022).
Soil erosion by water is often initiated by raindrop splashes striking surface aggregates where there is no protective vegetative cover, breaking aggregates down, and splashing soil particles in all directions, with geometry favoring a net downslope displacement (Ellison, 1947). The finer displaced sediments result in reduced surface infiltration capacity (Mohammed and Kohl, 1987; Ma et al., 2022; Rawls et al., 1990), and when precipitation intensity or snow melt rates exceed infiltration capacity, shallow overland flow mobilizes loosened surface particles in the runoff. Raindrop splash and shallow overland flow erosion are frequently combined and referred to as “interrill erosion.” When the overland flow reaches a certain rate, it merges into channels or rills, and the concentrated flow transports particles detached by splash and shallow flow as well as detaching particles from the rill bottom and sides due to the hydraulic shear forces or energy dissipation of the rill flow (Elliot, 1988). Particles in the rill flow that collide with the sides and bed of the rill detach additional aggregates and particles. This process of rill channel bed and side erosion is frequently referred to as “rill erosion.” When rill soil displacement is such that a layer of less erodible soil is reached, rills tend to widen as the channel sides erode and form what is often referred to as an “ephemeral gully” because the channel will likely be obliterated by subsequent tillage operations (Douglas-Mankin et al., 2020).
The Water Erosion Prediction Project (WEPP) Model
The WEPP Model is a physically based runoff and erosion prediction model developed by the USDA Agricultural Research Service (ARS) along with other federal agencies and universities. Early development was led and coordinated by an interagency Core Team (Foster and Lane, 1987). WEPP generally runs on a daily time step with average annual erosion predicted after simulations of 30 to 100 years with stochastically generated climate inputs. The WEPP model simulates crop growth, including biomass accumulation, canopy development, plant senescence, planting, harvesting, and residue accumulation and decomposition (Flanagan and Nearing, 1995). Additionally, the model computes a daily water balance that includes precipitation, soil evaporation, plant transpiration, infiltration, percolation, runoff, and subsurface lateral flow (Savabi and Williams, 1995). When rainfall intensity exceeds infiltration capacity, WEPP estimates runoff, water erosion, and sediment delivery. The WEPP model estimates erosion due to interrill and rill sediment detachment, transport, and deposition processes and, in watershed mode, channel sediment detachment, transport, and deposition (Flanagan and Nearing, 1995). In the 1989 version of WEPP, the equation for estimating interrill erosion was (Laflen et al., 1991; Nearing et al., 1989):
(1)
where
D_i = interrill detachment rate, kg m^-2 s^-1
V = vegetation factors
S_f = slope factor
K_i1 = interrill erodibility for equation 1, kg s m^-4
i = rainfall intensity, m s^-1.
In about 1989, WEPP team leaders implemented a change to the interrill detachment equation to better account for the combined effects of rainfall intensity and shallow overland flow (Foster et al., 1995).
(2)
where
K_i2 = interrill erodibility for equation 2, kg s m^-4
q = runoff rate, m s^-1.
Rill erodibility in WEPP is estimated from the shear of concentrated flow as described in detail in Elliot (1988), Elliot et al. (1989b), Flanagan and Nearing (1995), and Nearing et al. (1989). It varies spatially within a rill as hydraulic shear increases with distance downslope, but is expressed simply:
(3)
where
D_c = the rill detachment capacity for clear water at a given point in the rill, kg s^-1 m^-2
K_r = rill erodibility, s m^-1
t = hydraulic shear of the rill flow for a given point in an eroding rill, Pa
t_c = critical shear of the soil that must be exceeded by the rill hydraulic shear before rill detachment occurs, Pa.
However, as the sediment flux in the rill increases, the ability of the rill flow to detach sediment decreases, so the detachment rate D_r (kg s^-1 m^-2) describes the ability of sediment-laden rill flows to detach additional sediment, and in WEPP is modeled as:
(4)
where
G = sediment load at a given point in the rill, kg s^-1 m^-1
T_c = sediment transport capacity at the same point as G, kg s^-1 m^-1.
Transport capacity T_c is a function of hydraulic shear as described in Elliot (1988), Elliot et al. (1989b), Foster et al. (1995), and Nearing et al. (1989).
To develop the WEPP technology, databases for crop growth, tillage effects, conservation practices, and soil properties needed to be built (Foster and Lane, 1987). The widely used Universal Soil Loss Equation (USLE; Wischmeier and Smith, 1978) had a single value for soil erodibility, the USLE K Factor. WEPP, however, required three parameters for soil erodibility, plus additional soil properties for up to ten layers in the soil profile. The three erodibility parameters were baseline interrill erodibility (K_i) (eqs. 1 and 2), baseline rill erodibility (K_r) (eq. 3), and baseline critical shear (t_c) (eq. 3). The main WEPP hydrologic input parameter is baseline effective hydraulic conductivity (K_b; Alberts et al., 1995; Flanagan and Livingston, 1995). Additional soil-related coefficients to estimate sediment transport capacity and other hydrologic processes were determined within the model (Flanagan and Nearing, 1995). To estimate the erodibility values for the diversity of soils within the U.S., the WEPP Core team initially identified 30 cropland soils on which to measure the erodibility parameters (Foster and Lane, 1987). The selected soils comprised six soil orders, wet and dry climates, and cool and warm climates (Foster and Lane, 1987). The set included soils that were formed by aeolian, glacial, alluvial, and climatic weathering processes (Boul et al., 2003). Selection of specific sites for the soils of interest was often in the vicinity of ARS laboratories (Iowa, Nebraska, Oklahoma, Texas, Washington, Idaho, Minnesota, Missouri, Mississippi, Georgia, Maryland, and Indiana). Where this was not possible, state NRCS offices were contacted for assistance in identifying suitable sites for soils with the desired properties (Wyoming, California, Montana, North and South Dakota, Maine, New York, and Ohio). Similar diverse sets of soils were identified for rangeland and forest erodibility studies (Foster and Lane, 1987). It was believed that from these soils, regression relationships between readily measured soil properties and the three erodibility parameter values could be derived. From these regression equations, the erodibility parameter values could be estimated for many of the 20,000 soil series in the United States. The set of cropland soils was later expanded to 36 soils (table 1, fig. 1). Figure 2 shows samples of 34 of the soils sorted by geographic location, where the Ultisols in the southeast U.S. are bright red, as is soil 8 in Oklahoma. The Vertisol in Texas (Number 9) is black, as is the silty clay loam soil in South Dakota (Number 17). The Mollisols that dominate the northern U.S. are also noted for their darkness (soils in the center of the top three rows). Soil 13 on the left side of the second row is very light, likely due to its high calcium carbonate content.
Table 1. Site number, soil series, location, erodibility, and textural properties of WEPP cropland soils. Definitions of variables and their units are presented in table 2.
No. Soil series Location Texture K_i1^[a] K_i2 K_r t_c Clay Silt Sand FSi VfSa MSa Rock OC
1 Clarion Ames, IA FSL 2.03 2.03 4.6 0.4 18.9 27.2 53.9 14.4 10.9 13.8 4 1.05
2 Monona Castana, IA SiL 1.78 1.78 7.6 2.8 20.1 74.8 5.1 27 4.9 0 0 1.02
3 Cecil 1^[b] Watkinsville, GA SaL 1.18 1.18 8.4 2.2 8.6 16.7 74.7 11.5 5.9 22.4 4 0.55
4 Sharpsburg Lincoln, NE SiCL 1.65 1.85 5.3 3.2 39.8 55.4 4.8 23.9 4.6 0 0 1.85
5 Hersh Ord, NE SaL 3.83 3.93 11.2 1.7 9.6 13.4 77 2.2 32.9 18.2 0 0.49
6 Keith Albin, WY SiL 3.37 19.3 31.8 48.9 11.3 44 1 0 0.91
7 Amarillo Big Spring, TX SaL 3.82 4.12 45.3 1.7 7.3 7.7 85 2.9 21.1 23.9 0 0.16
8 Woodward Woodward, OK L 3.88 4 25 1.3 12.3 39.9 47.8 9.1 39 1.1 0.3 0.82
9 Heiden Waco, TX C 2 1.7 8.9 2.9 53.1 38.3 8.6 29.3 4.5 1 0 1.36
10 Whitney Fresno, CA SaL 2.66 2.74 23.3 4.7 7.2 21.7 71.1 9.7 8.1 19.2 4 0.19
11 Academy Fresno, CA SaL 3.06 2.88 5.7 1.6 8.2 29.1 62.7 12.2 20.2 13.7 4 0.41
12 Los Banos Los Banos, CA SiL 2.53 2.5 0.6 2.1 43 41.0 16 21.1 11.3 0.8 0 1.47
13 Portneuf Twin Falls, ID SiL 1.26 1.26 10.6 3.1 11.1 67.4 21.5 29.4 19.3 0.5 0 0.72
14 Nansene Dusty, WA SiL 3.13 3.12 30.7 3.1 11.1 68.8 20.1 30.6 18.1 0.3 0 1.49
15 Palouse Pullman, WA SiL 3.74 4.32 6.6 0.7 20.1 70.1 9.8 30.9 8.8 0.1 0 1.76
16 Zahl Bainville, MT L 3.21 3.17 12.3 3.5 24 29.7 46.3 14.7 12.5 7.2 9 1.69
17 Pierre Wall, SD SiC 2.17 2.18 11.7 4.8 49.5 40.9 9.6 21.4 7.3 0.8 0 1.46
18 Williams McClusky, ND L 2.95 2.94 4.5 3.4 26 32.4 41.6 17.5 11.5 8.7 5 1.79
19 Barnes Goodrich, ND L 1.72 1.71 3.3 2.3 24.6 36.0 39.5 20.6 12.1 6.7 4 3.26
20 Sverdrup Wall Lake, MN SaL 2.14 2.11 10 1.4 7.9 16.8 75.3 9.4 3.7 27.3 0.3 1.28
21 Barnes Morris, MN L 1.6 1.6 6.3 4 17 34.4 48.6 16.5 11.4 10.2 6 1.98
22 Mexico Columbia, MO SiL 2.97 2.97 3.6 0.7 26 68.7 5.3 33.1 1.1 1.5 0 1.56
23 Grenada Como, MS SiL 2.64 2.63 7.3 4.5 20.2 77.8 2 36.4 1.5 0.2 0 1.27
24 Tifton Tifton, GA Sa 0.76 0.77 11.3 3.5 2.8 10.8 86.4 4.8 13.3 26.2 23 0.46
25 Bonifay Tifton, GA Sa 1.58 0.87 17.9 1 3.3 5.5 91.2 3.1 16.2 18.6 1 0.32
26 Cecil 2^[b] Watkinsville, GA SaL 1.92 1.86 3.8 4.5 19.8 15.6 64.6 9.3 5.9 18.5 6 0.70
27 Hiwassee Watkinsville, GA SaL 1.88 1.88 10.3 2.3 14.7 21.6 63.7 15.3 4.3 19.2 3 0.83
28 Gaston Salisbury, NC CL 2.19 2.04 4.9 4.4 39.1 25.4 35.5 17.6 7.5 10.1 0.3 1.12
29 Opequon Flintstone, MD CL 3.3 3.2 3.5 6.3 31.1 31.2 37.7 22.8 5.9 9.6 14 1.42
30 Frederick Hancock, MD C 2.43 2.48 8.4 6.6 16.6 58.3 25.1 39.4 5.2 6.6 14 1.32
31 Manor Ellicot City, MD L 2.68 2.69 5.4 3.6 25.7 30.7 43.6 23.4 7.1 10.6 8 0.96
32 Caribou Presque Isle, ME L 1.55 1.45 4.5 4.3 12.2 40.8 47 25.6 11.5 7.8 47 2.28
33 Collamer Ithaca, NY SiL 3.74 3.46 24.13 6.4 15 78.0 7 47.5 4.6 0.7 0.3 1.01
34 Miamian Dayton, OH L 1.59 1.65 9.6 5.5 25.3 44.1 30.6 28.6 6.4 7.4 3 1.75
35 Lewisburg Columbia, IN CL 2.42 2.47 5.9 3.4 29.3 32.2 38.5 19.4 10.9 7.2 6 0.87
36 Miami Waveland, IN SiL 2.03 1.97 9.5 3.3 23.1 72.7 4.2 36.9 2 0.8 0 0.82
^[a] Units for all properties are presented in table 2
^[b]“Cecil 1” was a “non-eroded” site, and “Cecil 2” was an “eroded” site, but both soils were classified as “Cecil”
The WEPP core team decided to employ rainfall and runoff simulations to measure interrill and rill erodibility and critical shear. The cropland erodibility field studies were carried out in 1986, 1987, and 1988, with the measured interrill, rill, and critical shear values published in Elliot et al. (1989b), Laflen et al. (1991), and Flanagan and Livingston (1995). Concurrent with the erodibility studies, baseline hydraulic conductivity was estimated from two small, covered plots in 1987 and four covered plots in 1988 on each rainfall simulation site (fig. 3; Elliot et al., 1989b; Laflen et al., 1991; Rawls et al., 1990).
Dr. John Laflen, a United States Department of Agriculture, Agricultural Research Service (USDA ARS) research engineer in Ames, Iowa, was the principal investigator on the cropland erodibility study. Laflen was appointed research leader of the USDA ARS National Soil Erosion Research Laboratory in West Lafayette, IN, in 1987, in the middle of the field study. With his new responsibilities, he was unable to find the time to carry out the in-depth analyses presented by his team (Elliot et al., 1988, 1989a, 1990a, 1990b, 1993) or publish the relationships between soil erodibility and soil properties before he retired in 2000. Laflen did publish the erodibility results along with field methods in Laflen et al. (1991) when the cover of that journal issue was similar to figure 3. The authors felt that this study and its findings were of sufficient scientific interest and deserved to be published as a peer-reviewed article to support the soil erosion research community worldwide.
Figure 1. Location of 36 sites for 1986-88 WEPP Cropland Erodibility study (Laflen et al., 1991). Numbers refer to site numbers in table 1.
Figure 2. Range of colors of WEPP Cropland soils. Numbers refer to soil numbers in table 1. Soil Numbers 11 and 29 are missing. Soils are sorted by approximate geographic location.
Soil Erodibility
Soil erodibility is a quantification of the rate at which soil can be detached and transported by detaching or transporting processes of water or wind (Ellison, 1947; Lal and Elliot, 1994). The literature on the effect of soil properties on erodibility is considerable, and only some of the key references are reviewed in this article. One of the earliest summaries was by Bennett (1939). Because of his interest in the weathered soils of the southeastern U.S., chemical properties associated with weathered soils such as Fe and Al and base chemical (Na,Ca, and Mg) contents were evaluated along with texture and organic matter content. Smerdon and Beasley (1961) evaluated 11 Missouri soils and found that t_c was highly correlated with plasticity index (PI) and dispersion ratio, both of which are associated with clay content and mineralogy (Lambe and Whitman, 1969). Barnett and Rogers (1966) carried out rainfall and runoff simulations on 50 sites with 17 different soil series in the southeastern U.S. and found that combinations of textural properties, soil depth, and soil water content at field capacity were important predictors of sediment delivery, as was the plot slope steepness. Wischmeier and Mannering (1969) published a 24-term equation for estimating soil erodibility based on observations from 55 soils. Their equation included textural terms and ratios, organic matter, an aggregation index, soil water content and pH, variation of properties within the soil profile, slope shape, and loess effects. Wischmeier and Smith (1978) simplified the 1969 equation by using only soil texture, organic matter, and categories describing structure and permeability to estimate soil erodibility as the USLE K-Factor. Barnett and Rogers (1966) and Wischmeier and Smith (1965, 1978) assumed that soil erodibility described not only soil particle detachment and transport, but also the ability of the soil to generate runoff. El-Swaify et al. (1982) reported that soils that were more weathered (Ultisols and Oxisols) tended to have lower erodibilities than less weathered Inceptisols, Mollisols, or Alfisols. Meyer and Harmon (1984) reported on a rainfall simulation study of the interrill erodibility of 18 southeastern and midwestern U.S. soils, where they found that textural, organic carbon and chemical properties were correlated with observed interrill erosion rates. Meyer and Harmon (1984) also observed that the pH of the eroding water was correlated with interrill erosion.
Figure 3. Aerial view of WEPP cropland field experiments on site 17, with six rill plots adjacent to simulator (4 on left and two on right) with water flowing upward, and eight interrill plots on outside with ditches to drain away runoff from interrill plots. Elliot is on right side of simulator collecting runoff samples during Period 2 (rain plus added flow) from metal trays over interrill plots measuring potential runoff differences due to rain shadow effects by plot borders. Runoff samples from interrill plots were collected during Period 1 (rain only). Photo by Tim McCabe, USDA-ARS.
Onstad and Young (1982) stated that textural properties and aggregate stability were important predictors of both interrill and rill erodibility. They concluded that properties that contributed to strong aggregates, such as organic matter, bulk density, dispersion ratio, and moisture tension, were also important predictors of soil erodibility. Le Bissonnais (1996) summarized studies that found that soils with more stable aggregates were generally less erodible. He noted that aggregate stability was associated with the clay content and mineralogy, organic matter, and water content prior to stability measurements. Soils with a higher clay content or containing clay particles that were less likely to expand when wetted had more stable aggregates. In drier regions, the presence of sodium decreased aggregate stability, often expressed as the sodium adsorption ratio (SAR) between sodium ions and calcium and magnesium ions (Rahimi et al., 2000). The rill-interrill water erosion model in WEPP was first coded in the CREAMS model (Knisel, 1980). The CREAMS model used the USLE K-Factor in both the interrill and rill detachment algorithms. It combined the K-Factor with rainfall, runoff, topographic, cover management, and conservation practice factors similar in form to the USLE, but with modifications to both the rill and interrill detachment equations to account for differences in soil detachment processes (Knisel, 1980).
In summary, there is a large repository of articles discussing the relationships between soil properties and erodibility, but prior to the 1986-88 WEPP Cropland and Rangeland field erodibility studies, there were no published methods to specifically estimate the rill and interrill erodibility values required to apply the WEPP Model nationwide or elsewhere. The purpose of this article is to summarize the field procedures and describe the analyses that were carried out following the nationwide WEPP 1986-88 cropland field studies intended to develop regression relationships between the field-measured soil erodibility values, K_i, K_r, and t_c, and readily measured soil properties.
Methods
In the summer of 1986, rainfall and runoff simulation research was conducted on two soils in Iowa (Sites 1 and 2 in table 1 and also in fig. 1) to evaluate the erodibility data collection methods as described in Elliot et al. (1989b) and Laflen et al. (1991). A rotating boom rainfall simulator was used, with the nozzle height varying from 2.5 to 3 m above the sloping plots (Swanson, 1965). Four rill plots and six interrill plots were installed on each site. Rill cross-sections were estimated with a rillmeter (Elliot et al., 1997). Later in 1986, another study was carried out in Georgia (Site 3, table 1 and fig. 1) to test the use of stereo photogrammetry to evaluate rill cross-sections (Elliot et al., 1997). The WEPP Core team decided for the subsequent studies to install six 9-m long rill plots, six 0.5 x 0.75-mridged interrill plots, and two 0.5 x 0.75-m flat interrill erosion plots with the same slope as the site with no cover. In addition, two 0.5 x 0.75-m flat infiltration plots with a burlap cover were installed to estimate saturated hydraulic conductivity for 1987, as described in Elliot et al. (1989b), Laflen et al. (1991), and Rawls et al. (1990) (fig. 3).
Table 2. Definitions of soil properties and their units for the 1986-88 WEPP cropland erodibility study that were selected by at least one of the stepwise procedures for inclusion in regression equations to estimate soil erodibility in Elliot et al. (1988, 1989a, 1990a, b, 1993).
Variable Definition Units
15Bar Water content of soil at 15 Bar suction pressure often considered a measure of wilting point soil water content for the fine fraction (<2 mm dia) Percent
AgStab Aggregate Stability for aggregates < 5 mm Percent
Al Aluminum content in fine fraction Percent
Ca Calcium content in fine fraction MEQ/100 g
CaCO₃ Calcium carbonate content in fine fraction Percent
CEC Cation exchange capacity MEQ/100 g
Clay Clay content (<0.002 mm dia) of fine fraction (<2 mm dia) Percent
Cond Electrical conductivity of fine fraction mmhos/cm
Depth Depth of A Horizon cm
Fe Iron content in fine fraction Percent
FSi Fine silt content (0.002 mm – 0.02 mm dia) of fine fraction (<2 mm dia) Percent
K or K-Factor Soil erodibility value used in the USLE and RUSLE models. English units used in analyses (ton acre h [hundreds of acre-ft tonf in.]^-1)
K_i1 Interrill erodibility prior to 1989 x 10^-6 kg s m^-4
K_i2 Interrill erodibility from 1989 to present x 10^-6 kg s m^-4
K_r Rill erodibility x 10³ s m^-1
LL Liquid limit water content for soil fraction < 0.4 mm Percent
M Wischmeier and Smith’s (1978) texture property = (Silt + VfSa) x (100-Clay) Percent²
MC Soil water content measured in the field prior to the start of rainfall simulation g water/g soil
Mg Magnesium content in fine fraction MEQ Mg/100 g
Miner Mineralogy of the dominant clay fraction: 1 – Kaolinitic; 2 – Mixed; 3 – Mixed-Smectitic; 4 – Smectitic and Montmorillinitic; 5 – Calcitic (Elliot et al., 1990b)
MSa Medium sand content (0.25 mm - 0.5 mm dia) of fine fraction (<2 mm dia) Percent
N Nitrogen content in fine fraction Percent
Na Sodium content in fine fraction MEQ Na/100 g
OC Organic Carbon in fine fraction (< 2 mm dia) Percent
OM Organic matter in fine fraction (< 2 mm dia) Percent
Order Taxonomic Order of soil: 1 – Vertisols and Aridisols; 2 – Ultisols; 3 – Mollisols, Inceptisols, Alfisols, Entisols and Spodosols (Elliot et al., 1990b)
PL Plastic limit water content for soil fraction < 0.4 mm Percent
PI Plasticity Index = LL - PL Percent
Rock Rock content of whole soil sample (> 2 mm dia) Percent
Sand Sand content (0.05 mm – 2 mm dia) of fine fraction (<2 mm dia) Percent
SAR Sodium Adsorption Ratio (Na/((Ca+Mg)/2)^0.5); Measured on soils 1, 2, 4-22, and 33. Otherwise set to default value of 0.1 -
Silt Silt content (0.002 mm – 0.05 mm dia) of fine fraction (<2 mm dia) Percent
Slope Steepness of the study site. Percent
SpSf Specific surface of soil particles measured with ethylene glycol monoethyl ether (EGME) mg g^-1 EGME
t_c Critical shear of soil Pa
VSJ Handheld vane shear measurement of soil strength on field plot Pa
VfSa Very fine sand content (0.05 mm – 0.01 mm dia) of fine fraction (<2 mm dia) Percent
WDClay Water dispersible clay in fine fraction Percent
WDSilt Water dispersible silt in fine fraction Percent
The field methods and data analysis were described in Elliot et al. (1989a, b) and are summarized here. At each site, the rainfall simulator with Veejet 80100 nozzles was set up and leveled at the center of the research area. Plot slopes were determined with a surveying level. Six furrows were dug on each side of the simulator with a small tractor, with a ridging tool oriented across the slope if necessary to result in rill slopes of 3% to 6% for six 9-m-long rill plots. The other six furrows were dug for interrill plots (fig. 3). Metal collectors and pipes at the downstream end of the rill plots routed rill runoff to manual collection points. Metal borders were installed around the 0.5 x 0.75 m interrill plots, which were sloped toward the center following ridging, and a metal trough was cut to match the insloping plots installed at the bottom of each interrill plot. The trough drained to a small pit on the downhill side of the plot where timed grab samples were collected in 1-L plastic bottles. The four infiltration plots were made by leveling the furrows, installing metal borders on three sides of the plot, and a gutter on the downhill side of the plot to direct runoff to 1-L plastic collection bottles. Two of the infiltration plots were bare, and two were covered with burlap in 1987. In 1988, the scientist analyzing the runoff data requested that all four flat plots be covered with furnace filter materials from a local hardware store as the insloped interrill erosion plots were found to respond similarly to the bare flat plots in 1987, so all four flat plots could be used for estimating infiltration rates under a surface cover. Each interrill and infiltration plot had a manual rain gauge installed. During Period 1, simulated rainfall at ~60 mm h^-1 was applied, and runoff samples from the interrill, infiltration and rill plots were collected every five minutes until flow equilibrium was measured at the outlet of the rill plots. The total time for Period 1 varied from 30 to 60 minutes for most sites. Stereo photographs were also taken when the runoff from the rills stopped. Rill cross sections were then measured with a rillmeter, and soil strength measurements were carried out on the rill sides and bottoms with a pocket penetrometer, a handheld torvane shear device, and on the rill bottom, a fall cone penetrometer. During Period 2, rainfall simulation was resumed until the equilibrium runoff rate observed at the end of Period 1 was achieved. Then, a timed runoff sample was collected to measure runoff and sediment delivery rates, and the velocity in two rills was measured by timing the leading edge of a fluorescein dye plume spiked with a saline solution to travel 6 m in two rills. The observed leading-edge velocities were reduced by about 20 percent to account for dispersion within the plume based on velocity measurements made with an experimental continuous flow salinity meter that determined when the peak salinity concentration occurred on one of the rills during the 1987 season. As soon as the velocity measurements were complete and stereo pair photos taken, a nominal 0.13 L s^-1 flow was added to the top of each rill. When the increased flow was observed by the samplers, two consecutive timed samples were collected, rill velocities were measured, and stereo-pair photos were taken. The same procedure was repeated for flow additions of 0.26, 0.39, 0.52, and 0.65 L s^-1. On several sites in 1987 during Period 2, runoff was also collected from metal trays with borders similar in height to the interrill plot borders to evaluate any rain shadow effects from border walls under the rotating boom simulator (fig. 3). When the last pair of samples for 0.65 L s^-1 were completed, rainfall and flow additions were halted. Rill cross sections were measured again with a rillmeter, a stereo pair of photos was taken, and soil strength was measured. For Period 3, the same procedure used in Period 2 was followed, but without any rainfall. A final set of rill cross sections, stereo photos, and soil strength measurements were completed. All runoff sample bottles were taken to the lab to be weighed, dried at 105° C, and reweighed to determine runoff and sediment delivery rates. The runoff and sediment delivery results from every rill and interrill sample bottle (more than 7,000 bottles) were archived for the public in Elliot et al. (1989b). A set of additional soil strength measurements were carried out on a plot outside of the rainfall simulation circle with a procedure that a field technician could carry out on soil with an unknown erodibility. An infiltration ring was forced into the tilled soil, and 20 L of water was poured into the ring. After the water had all infiltrated, soil strength measurements with the torvane device, the pocket penetrometer, and the fall cone penetrometer were carried out.
Interrill erodibility values were calculated by solving equation 1 for K_i1 from the final three or four samples that were collected from the interrill plots when runoff rates were generally stable, as reported by Elliot et al. (1989a) and Liebenow et al. (1990). In 1987, two of the infiltration plots were bare, and they were also included in the erosion analysis. The slope factor in equation 2 was applied to the data to remove the differences due to the plot slope (Liebenow et al., 1990). In 1990, equation 2 was solved for K_i2, and all interrill plot runoff observations were used to estimate K_i2.
Determining K_r and t_c was more challenging, and an abridged description of the procedure is presented herein. See Elliot (1988) or Elliot et al. (1989b) for a more detailed description. Only rill data from Period 2 (rain plus added flow) were used for this analysis. Equation 3 showed that the observed sediment delivery rate (D_r) was a function of the hydraulic shear in the rill t. Shear was calculated from:
(5)
where
? = specific weight of water (N m^-3)
r_h = hydraulic radius (m)
s = plot slope (m m^-1).
? was determined from the eroding water temperature that was measured at each site (Elliot et al., 1989b). To estimate r_h for each runoff sample, the runoff rate (m³ s^-1) was divided by the rill velocity (m s^-1) to obtain the cross-sectional area of flow in the rill (m²). An optimizing computer program was written to determine the depth of flow needed in the rill cross-sectional shape measured with the rillmeter to give that area, and from that depth, calculate the wetted perimeter. The cross-sectional area was then divided by the wetted perimeter to calculate the r_h value for each rill runoff rate. This step also determined the width of the rill for that flow rate (w_r[m]).
Equation 4 required the transport capacity T_c of the rill flow for each observation. For this analysis, a simplified version of the Yalin sediment transport equation was developed (Foster and Meyer, 1972; Elliot, 1988):
(6)
where B was a sediment transport coefficient. B was derived from the aggregate and particle size distribution of the eroded sediment using the full Yalin (1963) equation for each of the five size classes (primary sand, silt, clay, and small and large aggregates) (Elliot, 1988). The value of B for each soil was archived in Elliot et al. (1989b).
An inspection of equation 4 shows that the detachment rate, D_r is a function of the amount of sediment in transport, G, but mathematically, D_r can also be expressed as the change in G with distance down the rill (Elliot et al., 1989b):
(7)
Substituting equation 7 into equation 4 allows for incorporating the effects of sediment already in transport G(x) on the change of change in sediment in transport dG/dx at a given point in the rill:
(8)
The amount of sediment in transport at any point in the rill plot (G(x)) increases with distance down the rill plot, and so sediment detachment will decrease with distance down the plot. The value for T_c was assumed to be constant. Equation 8 described rill erosion only. If sediment from interrill erosion (E) was added to the equation, as was the case during Period 2 (rain plus inflow), equation 8 becomes:
(9)
The units on E were adjusted so that all the interrill erosion from the 0.5-m wide ridges was concentrated in the observed rill width w_r (m). Equation 9 is a differential equation that has no exact solution, so Elliot (1988) developed an approximate solution to apply to the 11 observations of runoff and sediment delivery from a given rill in Period 2. Elliot’s (1988) approximate solution, assuming constant shear along the length of the rill, was:
(10)
where w_r was the width of the rill, L was the rill length (9 m), Q_s was the observed sediment delivery rate from the rill plot (kg s^-1) for a given sample; D_c was calculated from equation 3 with initially assumed values for K_r and t_c. The D_r values for all the 11 rill sediment delivery amounts for a given plot were then used in a regression analysis with shear for each rill flow. The slope of that analysis was an estimate for K_r, and the interception of that regression line with the x-axis was t_c. A revised value for D_c was calculated from the newly estimated K_r and t_c values with equation 3, reinserted into equation 10, and the analysis was repeated until equilibrium values for K_r and t_c were determined. Three or four iterations resulted in equilibrium values for K_r and t_c with three significant digits. The iterative process was carried out using a script in a spreadsheet. In a few cases where equilibrium could not be achieved, D_c in equation 3 was estimated from:
(11)
K_r and t_c were then estimated with a simple regression of D_c vs shear. Individual estimates for K_i1, K_r, and t_c were averaged for each site and are summarized in table 1. K_r and t_c were also calculated for all soils during period 2 using the D_c-computed values from equation 11, and those results were also reported in Elliot et al. (1989b). Using equation 11 for determining K_r resulted in lower estimates for K_r as the effect of sediment in transport was not considered, especially for the most erodible soils. Estimates for K_i2 were carried out after Elliot et al. (1989b) in Laflen et al. (1991).
In 1987, soil erodibilities were measured on 18 sites in the western U.S. (Sites 4 through 21, table 1 and fig. 1). Elliot et al. (1988) presented the results of that study as a progress report. In the initial 1987 analysis for interrill erodibility, we assumed that the rainfall energy from the rainfall simulator was likely only 80% of the energy from natural rainfall at the same intensity (Meyer and McCune, 1958), so a 0.8 factor was applied to the intensity. After the initial analyses, the investigators determined that this adjustment was unnecessary, and all interrill data from 1987 were reanalyzed. Some of the units used for the initial 1987 analysis were not in S.I. units (e. g., rainfall intensity changed from mm min^-1 to m s^-1 and sediment mass changed from g to kg; Elliot et al., 1988), so all units were altered to S.I. units, and the 1986 and 1987 data were reanalyzed. The 1988 data were only analyzed in S.I. units (Elliot et al., 1989b). The results from the Elliot et al. (1988) presentation with only the 1987 sites are included in this article to demonstrate how differences in calculating erodibility and in soil databases can impact the soil properties that were identified as important for estimating soil erodibility. In 1988, erodibility was measured on the final 15 sites in the Eastern U.S. (Sites 22 through 36 in table 1 and fig. 1), and the full data set was available for developing the regression relationships between soil erodibility values and measurable soil properties.
At each site, the Soil Conservation Service (SCS, now the Natural Resources Conservation Service, NRCS) carried out a soil survey with a pit up to 2 m deep dug near the center of the site, and up to four additional shallower pits excavated near the edges of each site. Detailed laboratory analyses of each pit's physical, chemical, and mineralogical properties were carried out at the SCS National Soil Survey Laboratory (now the NRCS National Soil Survey Center) in Lincoln, NE (USDA-SCS, 1990; Appendix B). Additional analyses of the engineering properties were also carried out by SCS specialists to see if common soil strength measurement techniques were useful for estimating soil erodibility.
The measured soil erodibility values (K_i, K_r, and t_c) were input into spreadsheets along with respective soil properties to format them for statistical analyses. Soil number 6 was not used for rill erodibility analysis as the site had excessive vegetation on it when the research team arrived, and subsequent root wads in the rill plots severely restricted rill erosion rates. The site for soil number six was in Eastern Wyoming, USA (fig. 1), but was selected by a scientist in Lincoln, NE, USA, some 670 km to the east. The local farmer had agreed to leave the site as fallow in the Fall of 1986 but failed to ensure that winter vegetation regeneration was not excessive. Surface residue or vegetation on all the other sites was minimal, and tillage effects from the ridging tool used to establish the plots were sufficient to ensure that existing stubble, surface residue, or vegetation was not a factor in the study.
More than 60 different soil physical, chemical, and mineralogical properties from the SCS Soil Survey Laboratory were also entered into spreadsheets. The data were exported as ascii files and were input into SAS Software analyses (SAS Institute, 1982) on the Iowa State University mainframe computer. An initial analysis of the correlation between erodibility and soil physical and chemical properties was carried out by Elliot et al. (1989a), as well as a correlation analysis among the soil properties. A special note was made of strong correlations among soil properties, such as clay and cation exchange capacity (CEC). Some properties were combined to provide additional variables for consideration that may better describe soil, such as the CEC:Clay ratio as an indicator of clay mineralogy, resulting in more than 70 potential variables for regression analyses. A stepwise procedure was then carried out, and the SAS software selected the best three- to four-term regression equations to estimate soil erodibilities from soil properties. It was not unusual to have a third or fourth variable selected that was highly correlated with one of the first variables selected, but with an opposite positive or negative sign. In such cases, that second property was removed from the list of potential variables, and the stepwise procedure was rerun.
In Elliot et al. (1989a), the authors focused on relationships between the SCS soil properties and soil erodibility. In the Elliot et al. (1990a) analysis, the engineering properties from the SCS were added to the other USDA-SCS (1990) soil properties. The Elliot et al. (1990b) analysis added geomorphological properties to the analysis, like site slope steepness, taxonomic order, mineralogy, and climate attributes. The Elliot et al. (1993) paper attempted to develop a set of nomographs like the USLE K-Factor nomographs (Wischmeier and Smith, 1978) using a limited common set of soil properties to estimate the three erodibility values.
The WEPP leaders decided that the regression relationships presented in the Elliot et al. (1989a, 1990a, 1990b, 1993) proceedings papers were unnecessarily complicated, and so a much simpler but less accurate set of predictive equations was incorporated into the WEPP model as published by Alberts et al. (1995).
As part of the preparation of this manuscript, all predictive equations were reevaluated with correlation and regression analyses in Excel to recalculate the Correlation Coefficients (r) between observed erodibility values and soil properties, and the Coefficient of Determination (r²) for each regression equation. The Elliot et al. (1988) regressions were also redone in Excel for this article, retaining the same properties identified in the 1988 analyses as important, but using the erodibility values that were common in the analyses completed by Elliot et al. (1989a, 1990a, 1990b, 1993).
A new set of predictive equations was developed for this article using a limited number of variables (soil texture, organic carbon, slope steepness, and taxonomic order) for incorporation into the WEPP model. The simplicity of the equations now in WEPP (Alberts et al., 1995) was enhanced using some of the relationships presented in the Elliot et al. (1989a, 1990a, 1990b, 1993) equations to improve the goodness-of-fit. The limited number of variables allowed us to develop the new set of equations using the multiple regression functions within Excel.
Results
This project aimed to develop regression relationships to estimate WEPP soil erodibility values from measurable soil properties. More than 60 soil or related properties were considered for the analysis from the USDA-SCS (1990) soil survey observations. Only 32 properties or ratios were selected in the stepwise procedures for inclusion in at least one of the regression analyses presented in the four proceedings papers (table 2; Appendix A). The selected soil physical properties are presented in tables 1 and 3a, and the selected soil chemical properties are presented in table 3b. Properties that were not selected included average annual temperatures and precipitation depths, all the NRCS engineering tests for soil strength, including a pinhole test intended to predict the likelihood of piping within earthen reservoir embankments, and all of the field strength tests except the vane shear test (VSJ) carried out on a plot outside of the rainfall simulation circle.
In table 1, the rock content was zero on soils with aeolian (Nos. 2, 4, 7, 14, 15, 22, and 23) or ancient lakebed (9, 12, and 17) geologic histories. Glacially derived soils (Nos. 1, 16, 18, 19, 21, and 32) had greater rock contents as did colluvium-derived soils (Nos. 10, 11, 29, and 30). Soil 9 had the greatest clay content (53%), a lakebed soil in eastern Texas. Soil 33 had the greatest silt content (78%) in upstate New York, where the soil was derived from lakebed deposits eroded from what are now the Appalachian Mountains. The greatest sand content (86%) was found in soil 25 derived from marine sediments on the Georgia coastal plain. Soil 19 had the most organic carbon content (3.3%) found on a glacial till soil in North Dakota with one of the lowest average annual temperatures (5.2°C) in the data set, while soil 7 in Texas had the least organic carbon (0.16%) and was one of the warmer sites with an average annual temperature of 17.4°C.
In table 3a, the torvane shear measurements (“VSJ”) on the external plot were not available for the 1986 soil numbers 1, 2, and 3, nor the soil water contents prior to the start of the rainfall simulations. The high value for M for soil 33 (M = 7021) reflects this soil’s claim to have the greatest USLE K-Factor in the U.S. Two other soils in the data set had even higher M values (Nos. 13 and 14), but likely had lower K Factors because they had greater hydraulic conductivities (Flanagan and Livingston, 1995). The greatest plastic limit (PL, 29%) was on soil 32, the soil with the greatest rock content (47%). The PL was measured on the soil particles less than 0.4 mm in diameter (table 2), and the large PL measurement on soil 32 may have been due to the high organic carbon content of that soil (2.3%; Gui et al., 2021).
The soil with the greatest aggregate stability (61%) was number 9, a Vertisol in Eastern Texas. Soil 9 was also one of the driest soils at the start of the rainfall simulation (0.003 g g^-1, table 3a). During the field study on soil 9, the only Vertisol in the database, rounded soil aggregates 1 to 3 mm in diameter were observed rolling down the rills, as there were few finer particles clouding the water near the end of the study. The highly stable aggregates, however, did not appear to have resulted in particularly lower erodibility values (table 1), with a K_r value of 0.0089 s m^-1 for soil 9 compared to the study’s mean K_r value of 0.0104 s m^-1.
In table 3b, the soil with the greatest iron and aluminum concentrations was Soil 28, a highly weathered soil derived from colluvium in the Appalachian Piedmont region of South Carolina. The greatest sodium concentration was on soil 17, a lakebed silty clay in South Dakota. The SAR value (Sodium Adsorption Ratio) was only calculated for western soils and soil 33, a lakebed-derived soil in the eastern U.S. Where SAR values were not available, a default value of 0.1 was assumed for regression analysis to allow consideration of logarithms of properties or ratios in regression equations (table 3b). The greatest SAR value was 0.3 for 17 soils. The calcium carbonate concentration was zero on many soils but very high on soil 9 (45%), a marine clay that was derived from calcareous mudstone. Soil 13, the only aridosol in the data set, also had a high concentration of calcium carbonate (23%). The 9% CaCO₃ concentration measured on soil 35 appears to be an outlier compared to the surrounding soils, which are all glacial till soils in Indiana or Ohio.
Table 3a. Physical properties that were useful for predicting WEPP cropland soil erodibility values. The definitions of the variables and their units are presented in table 2. The taxonomic orders that were split for estimating K_i2 in equation 4a are highlighted.
No M WD Clay WDSilt PL LL PI 15Bar Ag Stab Miner Miner/Clay MC Slope Order VSJ
1 3090 3.1 38.1 18 29 11 8.8 7 3 0.159 5.17 Hapludoll
2 6368 2.9 58.4 23 37 14 10.2 9 3 0.149 12.62 Hapludoll
3 2066 7.5 19.3 15 18 3 3.6 10 1 0.116 4 Kanhapludult
4 3612 7.5 62.7 25 19 11 17.4 4 4 0.101 0.06 5.68 Argiudoll 51972
5 4186 1.1 18.4 17 37 14 4.8 9 3 0.313 0.067 6.61 Ustorthent 7844
6 6117 3.4 40.7 19 18 3 10.3 6 3 0.155 0.03 5.08 Argiustoll 27456
7 2670 3.5 10.6 17 40 15 3.4 3 2 0.274 0.11 3.59 Paleustalf 13238
8 6920 3.2 43 24 18 1 5.9 3 3 0.244 0.247 7.14 Haplustept 21573
9 2007 17.6 59.6 20 29 10 18.9 61 4 0.075 0.003 3.86 Chromustert 19121
10 2765 4.1 23.4 13 18 1 2.6 13 2 0.278 0 7.36 Haploxeralf 52167
11 4526 2.6 34.2 17 28 4 3.5 9 2 0.244 0.033 4.51 Haploxeralf 31379
12 2981 12.3 54.5 21 55 35 19.1 15 4 0.093 0.03 4.03 Haploxeralf 68543
13 7708 3.5 77.6 23 15 2 8.3 33 2 0.180 0.053 5.57 Haplocalcid 67661
14 7725 2.2 71.7 24 18 1 7.3 20 2 0.180 0.097 6.1 Haploxeroll 65209
15 6304 2 69.8 23 49 28 9 20 2 0.100 0.13 6.47 Haploxeroll 38243
16 3207 6.2 36.9 18 27 4 10.2 12 3 0.125 0.133 7.56 Calciustoll 56875
17 2434 14.3 53.3 14 28 4 19.4 40 4 0.081 0.18 6.65 Haplustept 45304
18 3249 8.8 44.8 19 32 9 11.9 5 3 0.115 0.15 5.05 Argiustoll 73545
19 3627 7.1 43.8 27 34 16 12.7 13 3 0.122 0.09 5.78 Hapludoll 55894
20 1888 2.8 21.7 23 24 10 5.4 13 2 0.253 0.093 4.23 Hapludoll 57855
21 3801 8.6 41.5 18 36 17 9.3 6 3 0.176 0.17 8.27 Hapludoll 55894
22 5165 6.4 82.8 23 40 13 11.5 6 3 0.115 0.1 3.87 Epiaqualf 17651
23 6328 5 89.9 23 33 10 10.9 5 3 0.149 0.01 8.68 Fraglossudalf 17651
24 2343 2.5 11.8 17 34 16 2.1 3 1 0.357, 0.065 4.53 Kandiudult 33340
25 2098 1.5 7.6 17 40 17 1.1 3 1 0.303 0.03 3.9 Paleudult 37262
26 1724 14.4 18.4 17 33 10 8.3 15 1 0.051 0.01 4.32 Kanhapludult 56875
27 2209 10.6 26.4 16 18 1 6.4 3 1 0.068 0.11 3.96 Kanhapludult 33340
28 2004 24.5 32 20 18 1 14.8 23 1 0.026 0.115 6.38 Hapludult 39224
29 2556 24.4 58 26 29 12 14.2 29 2 0.064 0.13 12.05 Hapludalf 29418
30 5296 8.4 60 26 25 9 6.8 18 2 0.120 0.115 12.85 Paleudult 37263
31 2809 12.6 41.6 23 39 19 11.7 6 2 0.078 0.065 8.54 Dystrudept 3137
32 4592 7.7 41.9 29 38 12 8.9 10 2 0.164 0.135 8.85 Haplorthod 23534
33 7021 6.9 81.7 23 30 8 6.9 10 2 0.133 0.025 8.68 Hapludalf 23534
34 3772 12.2 52.9 21 44 21 10.9 10 2 0.079 0.11 8.92 Hapludalf 31379
35 3047 19.4 40.6 16 32 3 11.3 18 1 0.034 0.025 7.45 Hapludalf 15689
36 5744 7.5 86.6 21 29 6 10 5 3 0.130 6.5 Hapludalf 23534
Table 4 lists some of the properties that were evaluated but not selected by the stepwise procedure for estimating any of the erodibility parameter values. However, some of the properties were correlated with soil erodibility, as noted in table 4. The temperature and sodium content of the eroding water had no effect on soil erodibility, according to the 1988 study. The Elliot et al. (1989a) study found that the coefficient of linear expansion, an indication of the potential swelling of the clay component of the soil and aggregate stability, was not useful for predicting erodibility values, nor was the water content at field capacity (1/3 bar). None of the laboratory-measured engineering properties of the soils evaluated in the Elliot et al. (1990a) study were useful for estimating soil erodibility (table 4). The fall cone penetrometer field measurement had the greatest correlation with K_i2 (r = -0.52) of the soil strength properties that were not selected in the stepwise procedures for inclusion in multiple regression.
Table 5 presents the correlation coefficients (r) between the 34 soil properties and the three ratios that were selected during the stepwise procedures with the soil erodibility values and the clay content for at least one of the Elliot et al. presentations. The clay content was correlated with many other soil properties and tended to confound the stepwise procedure that assumed that all the potential variables were independent. The correlation analysis revealed that K_i1 and K_i2 were highly correlated (r = 0.98), suggesting that the predictive regression equations developed for K_i1 in Elliot et al. (1988, 1989a) could be applied to K_i2 with little loss of accuracy. Other correlation coefficients in table 5 for K_i1 and K_i2 were similar for all the properties listed. Table 5 shows that the very fine sand (VfSa) content was the only property with a correlation coefficient greater than 0.25 or less than -0.25 for all three erodibility properties and clay content. K_i1 and K_i2 were correlated with K_r, M, rock content, and soil taxonomic order. K_r was positively correlated with K_i1, K_i2, and the mineralogy:clay ratio, and negatively correlated with clay and rock content, organic carbon, iron, aluminum, calcium, magnesium, specific surface area, and 15-bar (wilting point) water content. All the negative correlations for K_r were positively correlated with clay content, suggesting that some care was needed when evaluating K_r stepwise results. Critical shear was positively correlated with fine silt, rock, water-dispersible clay, iron, and aluminum contents, aggregate stability, and plot slope steepness, and negatively correlated with the Miner:clay ratio, Sand, and very fine sand contents (VfSa).
The variables that were selected for the stepwise regression analyses are presented in table 6 for each of the Elliot et al. (1988, 1989a, 1990a, 1990b, and 1993) proceedings papers and Alberts et al. (1995). The full equations are given in Appendix A. The properties listed in table 6 were often combined as ratios, products, the argument for the natural logarithmic function, or the exponent of the natural logarithm constant e in the final equations (Appendix A). For some of the predictive equations, the data were divided into two sets before carrying out the stepwise procedure. The most common variables selected were clay (9 equations), sand (8 equations), and VfSa (8 equations). WDClay was in 6 equations, often in a ratio with Clay. The equations with the lowest Coefficients of Determination were the Elliot et al. (1988) set, where only the soil properties and erodibility values from the 1987 field season were considered, and the Elliot et al. (1993) and Alberts et al. (1995) sets with limited properties to simplify the predictive equations (table 6).
Table 3b. Chemical properties for soils presented in table 1 that were useful for predicting soil erodibility in the WEPP cropland soil erodibility study. Definitions and units for the variables are presented in table 2.
No Fe Al N Ca Na Mg CEC SAR Cond CaCO₃ Sp Sf
1 0.9 0.1 0.106 12.1 0.1 3.6 14.1 0.3 0.48 0.0 19
2 1.2 0.1 0.112 17.2 0.1 4.8 18.9 0.1 0.46 0.0 19
3 0.9 0.1 0.049 2.3 0.03 0.7 3 0.1 0.6 0.0 4
4 1.1 0.2 0.171 19.4 0.1 6.6 29.4 0.3 0.2 0.5 36
5 0.2 0.1 0.057 5.6 0 1.5 7.6 0.3 0.04 0.5 7
6 0.5 0.1 0.093 12.9 0.03 3.1 18.3 0.1 0.38 0.5 15
7 0.2 0.03 0.045 5.9 0.3 0.5 5.1 0.3 0.25 0.5 6
8 0.6 0.03 0.092 0 0.03 1.9 10.4 0.3 0.56 1.0 8
9 0.6 0.1 0.125 0 0.1 1 33.3 0.3 0.66 45. 38
10 0.5 0.03 0.024 2.8 0.1 0.7 3.4 1 0.11 0.5 3
11 0.9 0.1 0.034 3.2 0 1.4 5.3 0.3 0.84 0.5 4
12 0.8 0.1 0.146 28.3 0.5 8.6 38.7 1 0.89 1.0 45
13 0.3 0.03 0.083 0 0.2 6.3 12.6 1 0.76 23. 11
14 0.9 0.1 0.112 8.6 0.1 2.3 15.1 0.3 0.93 0.5 11
15 0.9 0.1 0.149 12.7 0.1 2.5 19.6 0.3 1.37 0.5 18
16 0.9 0.1 0.173 22.6 0.1 4.3 20 0.3 0.72 0.5 23
17 1.5 0.1 0.16 32.5 0.1 7.1 36.8 0.3 0.5 0.5 46
18 0.7 0.1 0.173 24.6 0.1 6.2 22.7 0.3 0.81 0.5 25
19 0.9 0.1 0.298 13.8 0.03 5 24.1 0.3 0.42 0.5 21
20 0.5 0.1 0.108 9.5 0 2.1 11 0.3 0.35 0.5 8
21 0.8 0.1 0.181 0 0 5.1 19.5 0.3 0.72 1.0 15
22 1.6 0.2 0.152 15.6 0.1 3 21.3 0.3 0.5 0.5 34
23 1.7 0.2 0.129 4.1 0.1 1.5 12 0.1 0.3 0.5 17
24 0.5 0.1 0.034 1.5 0.03 0.3 2.1 0.1 0.2 0.5 3
25 0.2 0.1 0.015 0.7 0.1 0.1 1.7 0.1 0.2 0.5 3
26 2.4 0.2 0.057 2.2 0.03 0.6 3.5 0.1 0.2 0.5 8
27 1.5 0.2 0.075 1.6 0.03 0.4 4.4 0.1 0.2 0.5 7
28 4.5 0.5 0.107 4.8 0.03 1.8 9.2 0.1 0.2 0.5 18
29 2.7 0.3 0.16 10.3 0.1 1 12.9 0.1 0.2 0.5 18
30 1.4 0.2 0.135 5 0 1.1 8.4 0.1 0.2 0.5 10
31 3.1 0.4 0.108 7.9 0.03 1.7 11.9 0.1 0.2 0.5 17
32 1.3 0.4 0.193 5.5 0.1 1.5 12 0.1 0.2 0.5 15
33 1.1 0.2 0.105 7.8 0.1 0.5 8.9 0.3 0.2 0.5 10
34 1.4 0.1 0.153 10 0 4.5 12.5 0.1 0.2 0.5 21
35 1.6 0.2 0.099 11.2 0.03 3.2 14.4 0.1 0.2 9.0 17
36 1.4 0.2 0.091 11 0.1 1.2 13.3 0.1 0.2 0.5 20
Regression equations requiring a reduced set of input data were derived in support of this article. The regression analyses were limited to soil texture, organic matter, CEC, site slope, and taxonomic order; all variables that are currently required elsewhere in the WEPP model soil input file except taxonomic order. The following three predictive equations were derived using the multiple regression capabilities of Excel:
For K_i2, for Ultisols and Aridisols:
(12a)
For K_i2for all other Taxonomic Orders:
(12b)
(13)
(14)
where the variables and their units are described in table 2. The probability that a coefficient was zero was less than 0.05 for all coefficients, so all terms in equations 12a and 12b, 13, and 14 are statistically significant except for the intercept (0.266) in equation 12a (P(0.266=0) < 0.28). Figures 4, 5, and 6 show the erodibility values predicted by equations 12, 13, and 14 versus the observed values. The confidence limits for each regression equation are also shown (Snedecor and Cochran, 1972). These confidence limits, however, do not include the variability within the soil properties data set.
Table 4. WEPP cropland erodibility study soil properties evaluated but not selected in stepwise multiple regression procedures. Correlation coefficients (r) greater than 0.25 are noted.
Elliot et al. (1988)
Temperature on day of study (r = -0.31 with K_i2)
Final soil water content; Depth of A Horizon
Sodium content and density of rainfall and runoff simulation water
Elliot et al. (1989a)
Fine sand (r = 0.27 with K_r and -0.27 with t_c); Coarse sand (r = -0.31 with K_i2)
Very coarse sand (r = -0.32 with K_i2); USLE K-Factor (r = 0.357 with K_i2)
Soil water content at 1/3 bar (field capacity) (r = -0.32 with K_r and 0.27 with t_c)
Elliot et al. (1990a)
Coefficient of linear expansion (r = -0.30 with K_r)
Onsite torvane on rill bottom before rainfall (r = -0.40 with K_r)
Onsite torvane on rill side before rainfall (r = -0.41 with K_r)
Onsite torvane on interrill plot before rainfall (r = -0.29 with K_r)
Onsite torvane on rill bottom after rainfall (r = -0.34 with t_c)
Onsite torvane, pocket penetrometer and fall cone penetrometer in saturated soil outside of plots
Onsite pocket penetrometer reading of soil strength with big foot before rainfall (r = -0.26 with K_i2)
Onsite pocket penetrometer reading of soil strength with big foot after rainfall (r = -0.47 with K_i2)
Onsite fall cone penetrometer measurement of soil strength (r = -0.52 with K_i2)
Laboratory unconfined consolidation stress (r = 0.29 with t_c)
Laboratory soil stress at 1% consolidation (r = -0.36 with K_r)
and 4% consolidation (r = -0.33 with K_r and 0.34 with t_c)
Laboratory pinhole critical shear test (r = -0.36 with K_r); Laboratory angle of internal shear
Coefficient of Linear Expansion (r = -0.30 with K_r)
PI/Clay (r = -0.27 with K_i2, 0.36 with K_r and -0.31 with t_c); Sp Sf / Clay (r = -0.25 with t_c)
Clay – WDClay (r = -0.32 with K_r); (Clay – WDClay)/Clay (r = 0.42 with K_i2 and -0.42 with t_c)
Elliot et al. (1990b)
Average annual temperature of site; Average annual precipitation of site (r = -0.39 with K_i2)
Discussion
This article aims to describe the analyses that were carried out following the WEPP cropland erodibility field studies that were intended to develop regression relationships between the field-measured soil erodibility values, K_i, K_r, and t_c, and readily measured soil properties. The results showed that soil erodibility is correlated with soil textural, chemical, mineralogical, biological, and geological properties. These diverse properties and their interactions were utilized to predict erodibility parameter values.
Table 5. Correlation coefficients (r) between K_i1, K_i2, K_r, t_c, and clay content and soil properties in Elliot et al.'s (1989a, 1990a, 1990b, 1993) analyses for WEPP cropland soil erodibility studies. Soil property descriptions and units for variables are described in table 2. The soil properties with the greatest correlation coefficient for each erodibility value are highlighted.
Property K_i1 K_i2 K_r t_c Clay
K_i1 1.000 0.979 0.364 -0.171 -0.031
K_i2 0.979 1.000 0.405 -0.102 -0.034
K_r 0.364 0.405 1.000 0.049 -0.417
t_c -0.171 -0.102 0.049 1.000 0.196
Clay -0.031 -0.034 -0.417 0.196 1.000
Silt 0.125 0.166 -0.076 0.209 0.222
Sand -0.085 -0.118 0.252 -0.256 -0.634
FSi 0.034 0.063 -0.155 0.442 0.350
VfSa 0.415 0.368 0.249 -0.558 -0.327
MSa -0.210 -0.188 0.201 -0.070 -0.549
M 0.338 0.336 0.168 -0.083 -0.227
USLE-K Factor 0.337 0.357 0.183 -0.008 -0.166
Rock -0.284 -0.265 -0.189 0.314 -0.198
OC -0.121 -0.072 -0.408 0.156 0.437
OM -0.008 0.027 -0.395 0.186 0.410
WD Clay -0.117 -0.126 -0.350 0.495 0.698
WDSilt 0.151 0.174 -0.162 0.233 0.385
WDClay/Clay -0.406 -0.418 -0.002 0.451 -0.183
CEC/Clay 0.087 0.123 0.091 -0.360 -0.161
Fe -0.046 -0.032 -0.355 0.430 0.426
Al -0.082 -0.099 -0.373 0.373 0.291
N -0.042 0.009 -0.398 0.194 0.515
Ca 0.160 0.187 -0.256 -0.042 0.563
Na 0.128 0.209 0.196 -0.082 0.244
Mg -0.149 -0.121 -0.343 -0.061 0.509
CEC 0.019 0.025 -0.338 -0.073 0.808
SAR 0.073 0.127 0.150 -0.052 0.032
Cond 0.184 0.236 -0.075 -0.321 0.119
CaCO₃ -0.180 -0.205 -0.035 0.002 0.365
Ag Stab -0.070 -0.074 -0.099 0.256 0.579
SpSf -0.037 -0.028 -0.391 0.045 0.899
15Bar -0.051 -0.055 -0.477 0.174 0.964
PL 0.062 0.086 -0.171 0.177 0.101
LL 0.002 0.030 -0.174 -0.061 0.147
PI -0.033 0.026 -0.220 -0.065 0.147
MC 0.160 0.220 0.118 0.059 0.029
Slope 0.103 0.141 -0.096 0.601 0.058
Order 0.486 0.519 0.054 -0.088 0.025
Miner 0.191 0.211 -0.161 -0.151 0.559
Miner/Clay 0.071 0.053 0.490 -0.300 -0.731
VSJ -0.297 -0.275 -0.116 0.075 0.056
Six sets of regression equations were proposed and are presented in Appendix A. The set that was incorporated into the WEPP model (Alberts et al., 1995; Flanagan and Livingston, 1995) put a greater emphasis on using common properties in simple equations rather than achieving greater Coefficients of Determination by developing complex equations with less common soil properties. The WEPP Core team anticipated that from the many soil properties that were measured, a few soil properties would stand out as closely related to soil erodibility. The Correlation Coefficients in table 5 only had one such variable; very fine sand (VfSa) had correlations greater than 0.25 or less than -0.25 for all three WEPP erodibility parameters and was used in 8 of the 18 predictive equations. VfSa was first identified by Wischmeier and Mannering (1969) as a key property associated with soil erodibility. Itcontinues to be one of the key properties for estimating the RUSLE and RUSLE2 K-Factors (Renard et al., 1997; USDA-ARS, 2008; Wischmeier and Smith, 1978). The clay content was frequently selected (7 equations). Clay was often applied in ratios with water-dispersible clay or other clay properties that helped to describe the mineralogy of the clay (Appendix A). The importance of clay content was another factor identified by Wischmeier and Mannering (1969), as well as Barnett and Rogers (1966), who found numerous ratios of clay with other soil properties, that were useful for estimating soil loss. Like VfSa, Clay is included in the current estimates for the RUSLE and RUSLE2 K-Factors (Renard et al., 1997, USDA-ARS, 2008; Wischmeier and Smith, 1978).
Table 6. Soil properties selected by stepwise procedures and coefficients of determination (r²) for each stepwise multiple regression analysis in WEPP cropland soils erodibility study. Definitions and units for the variables are presented in table 2. Complete equations are in Appendix A.
Predictive variables to estimate soil erodibility
Analysis K_i1or K_i2 K_r t_c
Elliot al., 1988 for soil nos. 4-21 For K_i1: Clay, VSJ;
r² = 0.35 Clay, VSJ;
r² = 0.35 MC, Depth;
r² = 0.16
Elliot et al., 1989a For K_i1: Clay, AgStab, WDClay, Mg, Fe, Al, Cond; r²=0.73 M, CEC, SAR, Al, OC;
r²=0.71 Clay, MC, VfSa, CaCO₃, SAR, SpSf, WDClay, Sand; r²=0.61
Elliot et al., 1990a For K_i2: CaCO₃, Na, Ca, Mg, PI, Clay, SpSf, WDClay; r²=0.55 M, Sand, LL, Al, Mg;
r²=0.76 Cond, Sand, VSJ, WDClay, Clay, Silt; r²=0.59
Elliot et al., 1990b For K_i2: Miner, Cond, WDSilt, FSilt, CaMg, Na, Clay, WDClay, SpSf; r²=0.79 Miner, Clay, Al, Mg, 15Bar, M, VfSa; r²=0.76 Clay, WDClay, Sloper²=0.75
Elliot et al., 1993 For K_i2: Silt, VfSa, Sand, OC, PL, Mg; r²=0.35 Silt, VfSa, Sand, OC, PL, Mg; r²=0.81 Silt, VfSa, Sand, OC, PL, Mg;
r²=0.26
Alberts et al., 1995 For K_i2: Sand, VfSa, Clay;
r²=0.24 Sand, VfSa, OM;
r²=0.55 Sand, VfSa, Clay
r²=0.23
The role of soil chemistry in erodibility was clear with the predictive equations including Fe, Al, Na, Ca, Mg, and SAR (tables 5 and 6). Fe and Al were associated with more weathered soils in the southeastern U.S. (table 3b, Soils 26 through 31) and were negatively correlated with erodibility parameters. The cations, however, were more common in the younger soils in the Western U.S. Soil 13 was exceptionally high in CaCO₃, while soils 16 through 18 were high in Ca, and soils 12 through 18 were high in Mg. Ca and CaCO₃ that are associated with the cementing of soil aggregates (Buol et al., 2003), whereas Na and the sodium adsorption ratio (SAR) were positively correlated with K_i1, K_i2, and K_r, and are associated with more readily dispersible clay aggregates (Rahimi et al., 2000). Mg was negatively correlated with K_i1, K_i2, and K_r, suggesting that Mg may increase soil strength (Latifi et al., 2015), reducing erodibility. Electrical conductivity Cond was only measured on soils in the western U.S. and was used as an indicator of sodic soils that were easily dispersed. In this study, we found that Cond was positively correlated with K_i1 and K_i2, and negatively correlated with t_c, suggesting that Cond is associated with weaker aggregates, but Cond was not correlated with K_r. Mg was also positively correlated with Clay (r = 0.51), so its selection in the stepwise process may have been as a Clay surrogate.
Considering soil engineering properties in the Elliot et al. (1990a) analysis improved the K_rr² value, but not for K_i2 and t_c. The fall cone penetrometer had the greatest correlation coefficient of all the unused soil strength properties (table 4) with K_i2 (r = 0.52). Bradford et al. (1992) reported that fall cone-measured shear stress was correlated with splash erosion when combined with raindrop energy in a laboratory study. The fall cone penetrometer, however, did not easily lend itself to a field application, requiring the team to set up a stable platform to support the device on saturated and weak soil surfaces.
The Elliot et al. (1990b) properties set added plot steepness, soil taxonomic order, and clay mineralogy to the 1989a data set that was based on lab-measured properties only, increasing the r² values by 5 to 14%. The 1990b analysis had separate regression equations developed for K_i2 depending on clay mineralogy (eqs. A10a and A10b), and separate regression equations for K_r depending on the Miner/Clay ratio (eqs. A11a and A11b). This differentiation suggested that more weathered soils (Miner = 2) with kaolinitic clays were less erodible than younger soils associated with smectitic clays. For the smectitic soils (Miner = 3), Cond was an important predictor of K_i2, but not for the other soils. Cond was only measured on the smectitic western soils, so it was not a factor in the non-smectitic eastern soils. For the non-smectitic soils, anions, and the product of two clay ratios (Clay/SpSf and Clay/WDClay), were selected by the K_i2 stepwise procedure. The K_r predictive equations split by the Miner/Clay ratio were more complex, with a different set of physical and chemical properties important in predicting rill erodibility for the two subsets (eqs. A11a and A11b). It appears that soil chemistry and clay mineralogy are important soil properties associated with interrill erosion, whereas both soil physical and chemical properties along with mineralogy are associated with rill erodibility.
Figure 6. Estimated critical shear t_c from equation 14 and t_c confidence intervals vs observed critical shear for WEPP Cropland Erodibility study (Elliot et al., 1989b) with blue dots for estimated values, red dashes for upper confidence interval, and green dashes for lower confidence interval (a = 0.05). Predicted vs observed t_c value from Mirzaee and Ghorbani-Dashtaki (2021) is a yellow triangle. Solid line is 1:1 line. In legend, “e” is error estimate for confidence limits (Snedecor and Cochrane, 1972).
One of the unexpected findings was the high correlation between plot slope and critical shear (r = 0.6, table 5). This suggests that steeper slopes may have experienced greater erosion rates historically, resulting in the preferential removal of finer particles and leaving coarser, less erodible material on the slopes. There was some concern that the high correlation may be due to the analysis methodology because hydraulic shear was used to develop the estimate of critical shear and was calculated from the product of slope steepness and hydraulic radius (Elliot, 1988; Elliot et al., 1989b; Nearing et al., 1989). To address this concern, a study reported by Yao et al. (2008) evaluated critical shear on a sloping soil bed and found that for the same soil, the critical shear remained constant regardless of the bed steepness, confirming that critical shear was more likely a soil property based on the soil forming and geologic processes associated with slope steepness than an anomalous result from data analyses.
The simplified predictive set presented by Elliot et al. (1993), where the same six properties were used for all three erodibility parameters, worked well for K_r (r² = 0.81) but resulted in poor predictions for K_i2 and t_c. Elliot et al. (1993) may not have chosen the best six properties for their simplified predictive equations.
One soil property that was not selected very often in the analyses was the USLE K-Factor (tables 4 and 5). The English unit K-Factors in Laflen et al. (1991) were correlated with K_i2 (r = 0.36) but not with K_r(r = 0.18) nor t_c(r = -0.01). We calculated a significant correlation between the K-Factor and the baseline saturated hydraulic conductivity for the 36 soils presented in Flanagan and Livingston (1995) (r = -0.53), suggesting that the dominant process associated with sediment delivery from the historic USLE plots was not the detachability or transportability of sediment but rather the likelihood of a plot generating runoff. It also raises some concern about the validity of linking the K-Factor with a separate runoff estimate as done in CREAMS (Knisel, 1980) and USLE-based watershed models since the runoff is already accounted for in the K-Factor. The correlation between the USLE K-Factor and K_i2 (r = 0.36) suggests the dominant erosion process on the relatively short USLE plots from which the K-Factor was developed was interrill erosion.
Gilley et al. (1993) used the Elliot et al. (1989a) data set to develop predictive equations for critical shear. Their approach to estimating K_r and t_c, however, was not the same as that described by Elliot (1988), Elliot et al. (1989b), or Laflen et al. (1991), who plotted rill detachment vs shear to find the X-axis intercept for t_c. Gilley et al. (1993) plotted shear vs. detachment to find the Y-axis intercept for t_c. In the Gilley et al. (1993) regression analysis between t_c and soil properties, they grouped the soils by water dispersible clay (WDClay) and found that for soils with WDClay < 7.5%, critical shear could be predicted from the clay content, coefficient of linear extensibility (COLE), and soil water content at 1.5MPa (15 bar or wilting point); and for soils with WDClay = 7.5%, critical shear could be predicted from the potassium, Ca, Fe, OC, and VfSa contents. Three of the same variables (WDClay, VfSa, and Clay) along with Ca as part of the SAR, were selected in the Elliot et al. (1989a, 1990a, 1990b) equations, and VfSa and Clay were used in the Alberts et al. (1995) and Flanagan and Livingston (1995) regressions. In the Elliot et al. (1989a) analysis for t_c, they divided the soils by Clay content, but the properties they selected had more chemical properties in the lower clay content soils group (eq. A6b).
Bajracharya et al. (1992) reported an interrill erodibility study using rainfall simulation on five soils from Ohio. The K_i1 values they measured ranged from 1.23 to 1.97 x 10⁶ kg s m^-4 on low-gradient loam and silt loam soils. Bajracharya et al. (1992) used the Elliot et al. (1989a) K_i1 regression equation to estimate the interrill erodibility of these soils. A paired t-test on the Bajracharya et al. (1992) results indicated the predicted K_i1 values that were presented were not significantly different from the observed values (t = -0.17; P(T=t) = 0.87 for the two-tailed test). The predictive Elliot et al. (1989a) equation (eq. A4b) used some variables that were not readily available (Fe, Al, and Cond), but as there was little variation among these properties for soils in the region, default values were assumed from the nearby soils (soils 34, 35, and 36; table 3b). One of the soils in the Bajracharya et al. (1992) study was Miamian, the only Ohio soil in the study. For K_i1, Elliot et al. (1989b) measured 2.03 x 10⁶ kg s m^-4 compared to 1.97 x 10⁶ kg s m^-4 reported by Bajracharya et al. (1992), but only 0.86 x 10⁶ kg s m^-4 predicted by equation A4b for K_i1. The Miamian soil measured by Bajracharya et al. (1992) had a clay content of 18.1% compared to 25.3% in the WEPP Cropland data set. The lower Clay and likely lower WDClay contents (correlation coefficient is 0.698 between WDClay and Clay, table 5) in the Bajracharya et al. study were likely the reasons for the lower estimated K_i1 value. Table 5 shows that K_i1 was negatively correlated with Clay, WDClay, and especially the WDClay/Clay ratio (r = -0.41). The low predicted value for the Bajracharya et al. result also shows the sensitivity of equation A4b for estimating K_i1 to Clay and WDClay contents.
Soil 13 was problematic in predicting K_r because it had a high silt content but was also very high in carbonates (CaCO₃= 23%). It was the only aridic soil in the set. It was likely selected for the study because it is a widely distributed soil with high erosion rates under furrow irrigation in southern Idaho, U.S. (fig. 1; Bjorneberg, 2001; Buol et al., 2003). All the predictive equations except Elliot et al. (1993) overpredicted K_r by a factor of 2 to 3 for this soil, with the Alberts et al. (1995) equation overpredicting K_r by a factor of 5. Mirzaee and Ghorbana-Dashtake (2021) had a similar soil with high calcium content in Iran and found that the Alberts et al. (1995) equation overpredicted their observed K_r by a factor of 14. Soil 13 had the lowest K_i1 and K_i2 values of any of the silt loams and the greatest CaCO₃ concentration, suggesting that CaCO₃ may reduce interrill erodibility. Dimoyiannis et al. (2001) also observed lower interrill erosion in soils with greater CaCO₃ concentrations. When considering the entire cropland data set, however, the correlation coefficients for CaCO₃ with K_i2, K_r, and t_c were only -0.21, -0.04, and 0.002, respectively (table 5). It is likely that other soil properties tended to mask the effect of high CaCO₃ concentration values that were less than or equal to 1% on 33 soils, but were 45, 23, and 9% on soils 9, 13, and 35, respectively. This overprediction for soils high in CaCO₃ suggests that further research is warranted on such soils to better predict soil erodibility. Soils high in CaCO₃ are common in areas of low rainfall (Buol, 2003), and on coal surface mine reclamation sites (Eriksson and Daniels, 2021; Liu, 2012).
Although numerous references suggested that aggregate stability should be a good predictor of soil erodibility (Le Bissonnais, 1996; Wischmeier and Mannering, 1969; Onstad and Young, 1982), we found that the correlation coefficient values for K_i1, K_i2, K_r, and t_c were only -0.07, -0.07, -0.1, and 0.26, respectively. Only the critical shear showed some correlation with aggregate stability. The only regression equation with AgStab as a term was for estimating interrill erodibility in Elliot et al. (1989a) (eq. A4a) for soils with greater than 35% clay content. Meyer et al. (1992) observed in their study of interrill erodibility on 22 cropland soils in the southeastern and midwestern U.S. that eroded sediments consisted mainly of aggregates on all but the sandiest soils. Their results underscored the importance of aggregate stability in the erosion and sediment delivery processes, even if it does not appear to be a major factor in estimating soil erodibility. Meyer et al. (1992) also suggested that rill erosion and concentrated flow processes would likely be affected by the presence of aggregates. In the WEPP model, the lower density but larger diameter eroded aggregates make up two of the five sediment classes (Elliot, 1988; Flanagan and Nearing, 2000; Foster et al., 1995). The sediment transport capacity T_c in equation 3 incorporates the effects of larger diameter and lower density aggregates in its estimation (Elliot, 1988; Elliot et al., 1989b). Refer to Meyer et al. (1992) for an additional discussion on the detachment and transport of aggregates in water erosion.
For forest and rangeland conditions when applying WEPP (Elliot, 2004) or the rill/interrill Rangeland Hydrology and Erosion Model (RHEM) (Nearing et al., 2011; Al-Hamdan et al., 2017), the developers found that soil erodibility was influenced by the dominant perennial vegetation (forest, shrublands, bunch grasses, sod-forming grasses, and post-wildfire) and soil textural properties. Nearing et al. (2011) found that interrill erodibility depended on plant community (grasses and forbs or shrubs), vegetation canopy, litter, and rock cover amounts, but had insufficient data to estimate concentrated flow erosion. Al-Hamdan et al. (2017) found that for WEPP, K_i for rangelands could be estimated from the plant community (bunch grass, sod grass, or shrubs), litter, canopy and rock cover, soil texture, and effects of fire and tree encroachment. These results suggest that an additional factor to consider for improving cropland erodibility estimation might be pre-European land cover. The eastern soils generally had hardwood forests, the upper Midwest and Northwest soils, bunch grasses, and for the aridic soil and the California soils, shrubs. These pre-European cover effects, however, are likely reflected in the soil order variable, with forests associated with Alfisols, grasslands with Mollisols, and shrub and sparse grass rangelands with Aridisols (Buol et al., 2003). However, Alfisols and Mollisols had the same taxonomic order ranking in Elliot et al. (1990b; table 2), and efforts to improve the predictive performance of this article by adding additional taxonomic order categories or vegetation categories did not improve the regression equations for estimating soil erodibility.
The Alberts et al. (1995) predictive equations used soil textural properties plus organic matter (not organic carbon) in regression equations to estimate soil erodibility. They had the benefit of simplicity but the disadvantage of low r² values. The Elliot et al. (1989a and 1990b) sets of equations provided the best r² values, but some of the properties needed for the predictive equations are not readily available, like clay mineralogy, WDClay, Cond, 15Bar, or soil anion concentrations. When the study was initiated, the Core team directed that this study should look at as wide a range of soil properties as possible to find the best properties for predicting erodibility (Foster and Lane, 1987). As it turned out, simplicity was more important than accuracy in the final parameterization equations published as defaults in the WEPP model (Alberts et al., 1995; Flanagan and Livingston, 1995). The goodness-of-fit of the default WEPP equations could be significantly improved with the addition of taxonomic order for K_i2, reformulating the simple set of variables for K_r, and incorporating slope steepness for t_c (eqs. 12a & b, 13 and 14). The earlier studies (table 6 and Appendix A) suggest that equations 12a and 12b for K_i and equation 14 for t_c may be improved by incorporating terms with WDClay, and equation 13 for K_r may be improved by incorporating an Mg term.
The performances of the simplified equations 12, 13, and 14 were evaluated with independent data sets from Bajracharya et al. (1992) for K_i2 and Mirzaee and Ghorbani-Dashtaki (2021) for K_i2, K_r, and t_c. The results are plotted in figures 4, 5, and 6. All five measured Bajracharya et al. (1992) K_i2 values were less than the predicted values. The high estimated K_i2 values for the Bajracherya et al. study appear to be due to the relatively high VfSa content (17.6%) on one of the soils in that data set and a high silt content (61.5%) on another (fig. 4), increasing the estimated value for K_i2. Another possible reason for the overprediction of erodibility by equation 4b was that the rainfall simulator nozzle height used by Bajracherya et al. (1992) was only about 2 m, less than the 2.5 to 3 m needed for simulated rainfall to achieve terminal velocity (Meyer and McCune, 1958; Meyer and Harmon, 1979), perhaps leading to lower erosion rates for a given rainfall intensity. Both the K_i2 and the K_r values were overestimated for the Mirzaee and Ghorbani-Dashtaki (2021) soils. The observed t_c value determined by Mirzaee and Ghorbani-Dashtaki (2021) was 3.49 Pa compared to the estimated value of 3.48 Pa assuming their plots had average slopes of 5.56%, the average of the Elliot et al. (1989b) plots. Mirzaee and Ghorbani-Dashtaki attributed their lower erosion rates to the cementing effects of the CaCO₃ in the soil. The interrill rainfall simulator used in their study had a capillary-type device forming raindrops about 0.5 m above the plot, well below the drop height of the Elliot et al. (1989b) data set of 2.5 to 3 m (Swanson, 1965), so the raindrop energy was substantially less (Meyer and McCune, 1958; Meyer and Harmon, 1979), likely resulting in a lower measured interrill erosion rate. According to Elliot (2017), sediment concentration in runoff from rainfall simulators on interrill plots was 2.5 to 17 times higher from a nozzle simulator 2.5 m above the plots than from a drop-forming simulator 1 m above the plots. The rill erosion methods of Mirzaee and Ghorbani-Dashtaki (2021) measured the rill erosion rate with simulated runoff but no rainfall. The ARS WEPP cropland and rangeland erosion sets both had prewetting of the plots (Period 1), plus rainfall combined with runoff simulation when measuring rill erosion (Period 2). This resulted in soils being near saturation during the measurements of rill erosion for the WEPP datasets (Elliot et al., 1989b), but not for the Mirazaee and Ghorbani-Dashtatki data. Drier soils may have resulted in greater resistance to rill erosion due to soil water tension and hence lower K_r values (Elliot, 1988; Holz et al., 2015). Also, the lack of raindrop splash on shallow rill flow may also have reduced rill detachment rates (Kinnell, 1991). The Elliot et al. (1989b) data set included rill erosion observations for “flow only” (Period 3), but they did not calculate K_rand t_c values for the flow only data. It was apparent in Elliot et al. (1989b) that sediment delivery rates Q_s were lower from the “flow only” runs compared to the “rain plus flow” runs, although the reduced erosion may have been due in part to reduced availability of sediment in Period 3 as the flow only treatment was on the same rills that had been treated with the rain plus flow treatment, resulting in armoring on some soils (Elliot et al., 1989b; Foltz et al., 2008). The independent data sets underscore the challenges of estimating soil erodibility as soil properties can vary significantly within a given series, and different research methods and site conditions may result in different estimates for soil erodibility.
The rill erosion data were not fully analyzed for Period 3 in Elliot et al., 1989b. The observed reduced sediment delivery may have been due to increased soil strength due to increased soil water tension (Elliot, 1988; Holz et al., 2015), reduced soil availability due to armoring (Foltz et al., 2008), or a different rill morphology with lower flows winding their way down a widened rill channel from the previous period (Elliot, 1988). There is scope for additional studies to compare rill erodibility with and without added rainfall, as many recent studies use “flow only” to measure rill erodibility (Mirzaee and Ghorbani-Dashtaki, 2021; Pierson et al., 2008; Robichaud et al., 2010).
One aspect of soil erosion science that was not addressed in this study is variability. Elliot et al. (1989b) calculated the coefficient of variation (standard deviation divided by the mean) for the field-measured soil erodibility values that ranged from under 10% to greater than 40%, averaging 25%. With this level of uncertainty due to variation in a replicated study over a relatively small research site, trying to estimate soil erodibility values with a high degree of precision would be presumptuous. In addition to the variability in measuring soil erodibility, the variability of the measured soil properties has not been addressed. At each study site, the USDA-SCS (1990) collected five sets of soil samples: a central pit and four other pits nearer the perimeter of the 50 x 50 m study site. The variability of soil properties within each site has never been evaluated. Robichaud et al. (2007) addressed the variability in forest soils, ground cover, spatial variability, and climate. They incorporated those sources of variability into the post-wildfire Erosion Risk Management Tool (ERMiT) interface for WEPP for estimating the erosion associated with a given probability of occurrence (Al-Hamdam et al., 2022; Robichaud et al., 2007). Tiscareno-Lopez et al. (1995) described an approach for assessing uncertainty within the WEPP model for rangeland erosion prediction considering variability in soil texture, standing biomass, litter, and soil erodibility, and Brazier et al. (2000) did a similar analysis for WEPP on croplands considering the variability of WEPP inputs. No such analysis has been conducted for the errors associated with estimating cropland erodibility, considering the combined variability within the WEPP cropland erodibility study and the soil properties databases.
This study has identified some key variables for predicting soil erodibility, such as very fine sand, taxonomic order, organic carbon, and site steepness, and these variables are not likely to change if additional soils are added to the database. Additional analyses of the data set may reveal interactive properties or groupings among the soil properties that could improve the goodness-of-fit of the erodibility prediction equations. Adding soils may have the benefit of improving the predictive r² values if they do not differ widely from the soils already studied. However, the current Cropland Erodibility database has no Andisols, Histisols, Oxisols, or Gelisols, and only one Spodosol and one Aridisol (table 3a). Expanding and diversifying the database could increase the complexity of analysis, requiring more terms in predictive regression equations or more splits in the predictive equation sets to address increasing variability.
Conclusions
Rill and interrill erodibility values were measured on 36 cropland soils distributed across the USA. The goal of the study was to derive relationships to estimate interrill and rill erodibility and critical shear for cropland soils from measurable soil properties to be applied within the WEPP model. We presented six sets of equations that have been proposed to estimate interrill and rill erodibility and critical shear for cropland soils from measurable soil properties. In addition, we developed a reduced variable set of predictive equations. We investigated incorporating soil physical, chemical, engineering, mineralogical, topographic, climatic, and taxonomic properties into predictive regression equations. The very fine sand content was the only widely available property that was somewhat correlated to all three erodibility parameters. The WEPP model had, internally, the least complicated set of predictive equations with only soil texture and organic matter terms, but this set of equations gave the poorest estimates for rill and interrill erodibility coefficients of all the proposed predictors. The more complicated sets of equations that considered a wider range of soil properties performed better. We found that incorporating field- and laboratory-measured soil engineering properties did not improve the estimation of erodibility values. The best prediction equations included soil textural terms and terms that described the clay mineralogy or taxonomic classification. The critical shear term for estimating rill erosion was found to be highly correlated with the research site's steepness and a term describing clay content and mineralogy. We developed a set of erodibility estimation equations to consider incorporating into the WEPP model that included texture and organic carbon and a simplified method for accounting for clay mineralogy, the cation exchange capacity to clay ratio. We also recommend that the critical shear estimate include the slope of the site. Our findings suggest that further studies are required on soils high in calcium carbonate, as most of our erodibility estimation equations overpredicted the erodibility of the single silt loam soil in the data set that was high in calcium carbonate. There is also scope for interested researchers to use the extensive database of soil properties for these soils to develop alternative predictive equations that may improve on those that we have presented. There are research opportunities to carry out additional field studies to fill in gaps in the soil database or evaluate soils of local interest. Understanding the impacts of variabilities of soil erodibility and soil properties on erosion estimates also warrants further investigation.
Acknowledgments
We wish to acknowledge Dr. John Laflen (deceased) who was a USDA-ARS Research Agricultural Engineer when this study was initiated and was promoted to become the Director and Research Leader of the USDA-ARS National Soil Erosion Research Laboratory during the study. Dr. Laflen proposed the study described herein, supervised the collection of the field data and subsequent analysis, and made the final determination for estimating WEPP soil erodibility values. Dr. Laflen was also the major professor for both authors in their PhD studies. Dr. Kris Kohl deserves recognition for his role as a graduate student in collecting the cropland field data and in the subsequent analyses of those data. We acknowledge Mr. Richard Hartwig (ARS Research Technician), Ms. Ann Liebenow-Antilla (an undergraduate student who played major roles in data collection and analysis), as well as many other students who assisted with the field and laboratory data collection. Hartwig and Liebenow-Antilla were the only team members who were present at every site. We want to recognize Dr. Stephen Holzhey, former director of the SCS National Soil Survey Laboratory, Lincoln, NE, for providing the soil properties data, and the WEPP Core team for suggestions of research sites, research methods, and soil properties to consider. We thank David Hoover (Director) and Scarlett Murphy (Physical Science Technician) at the USDA-NRCS National Soil Survey Center, Lincoln, NE, for guidance on how to access the 35-year-old WEPP Cropland soil properties online database.
Appendix A
The regression equations derived for predicting WEPP cropland erodibility values were published as proceedings from five ASAE meetings (Elliot et al., 1988, 1989a, 1990a, 1990b, and 1993) and a sixth set in the WEPP User Summary (Flanagan and Livingston, 1995) and the WEPP Technical Documentation (Alberts et al., 1995). The equations are presented in this section. All variable definitions and units for these equations are presented in table 2.
Elliot et al. (1988) developed a preliminary set of equations as a progress report with data from soils 3 through 22 (table 1). The erodibility and property values they identified were updated following the presentation. Here are the three prediction equations using the original forms of the equations but with updated variables and recalculating the regression coefficients in Excel:
(A1)
(A2)
(A3)
The focus of the Elliot et al. (1989a) analysis was the large database of soil physical and chemical properties measured by the USDA-SCS National Soil Survey Laboratory (USDA-SCS, 1990).
For soils with Clay > 35%
(A4a)
For soils with Clay = 35%
(A4b)
(A5)
For Clay > 30%
(A6a)
For Clay = 30%
(A6b)
In 1990, the interrill erosion equation changed from equation 1 to equation 2, and K_i2 became the interrill erodibility parameter for all subsequent analyses. In Elliot et al. (1990a; eqs. A7, A8, and A9), the study focused on the field- and lab-measured soil strength properties but did not give the stepwise procedure the option of considering the geomorphic variables available in Elliot et al. (1990b). The 1990b analysis had the option of incorporating field-measured soil strength properties, however. The analysis was done in three ways: using field-measured strength values only, using lab-measured values only, and using both. The following equations resulted from the 1990a analysis that allowed the incorporation of both field- and lab-measured soil strength values. The result was that no lab-measured soil strength values were selected by the stepwise procedure (table 4). The equations presented in Elliot et al. (1990a) were:
(A7)
(A8)
(A9)
For the analysis reported in Elliot et al. (1990b), in addition to the USDA-SCS (1990) database, the authors considered soil properties related to taxonomic order, clay mineralogy from the SCS (1990) database, and the slope of the research site.
K_i2 for smectitic and mixed smectitic soils:
(A10a)
K_i2 for kaolinitic and mixed soils:
(A10b)
K_r for soils with Miner/Clay > 0.1 :
(A11a)
Kr for soils with Miner/Clay = 0.1 :
(A11b)
(A12)
The Elliot et al. (1993) equations were developed using a single set of properties for all three variables, and then developing a set of relationships that could provide a graphical method for estimating soil erodibility. Like the K-Factor nomographs for the USLE/RUSLE graphs (Wischmeier and Smith, 1978), there were four quadrants, with each quadrant incorporating another soil property (fig. A1). The following sets of equations have one equation for each quadrant. In the following equations, the variable SiFS was used to replace the sum of Silt + VfSa. Trigonometric function arguments are in radians.
(A13a)
(A13b)
(A13c)
(A13d)
(A14a)
Figure A1. Nomograph for estimating K_r x 10³ based on equations A14a to A14d. To use, start in upper left quadrant with appropriate textural amounts, then down to lower left quadrant to appropriate OC content, noting Krp1 on lower left axis. Move right from OC intercept, noting Krp2 to PL intersection. Move to bottom of lower right quadrant to note KrPp3, or up to upper right quadrant to intercept with Mg content, and then move right to estimate final value for Krp4, rill erodibility for soil of interest. Example shown in blue line is for soil 18 (Silt + VfSa = 45.8; Sand – VfSa = 37.2; OC = 1.98; PL = 18; and Mg = 5.1; Measured K_r= 6.3 x 10^-3 s m^-1 and estimated Krp4 = 7.3 x 10^-3 s m^-1).
(A14b)
(A14c)
(A14d)
When reevaluating the regression equations for t_c in Elliot et al. (1993), the r² values were very poor. The form of each equation was retained, but the coefficients for each of the terms were altered to achieve an equation set with better r² values for each equation. The nomograph was not redone.
(A15a)
(A15b)
(A15c)
(A15d)
The erodibility estimation equations in Alberts et al. (1995) and Flanagan and Livingston (1995),like those in Elliot et al. (1993), were based on a limited set of soil properties, including only soil textural fractions and organic matter (not organic carbon as used in previous equations).
For soils with Sand > 30%
(A16a)
For soils with Sand = 30%
(A16b)
For soils with Sand > 30%
(A17a)
where OM is the percent organic matter.
For soils with Sand = 30%
(A17b)
For soils with Sand > 30%
(A18a)
For soils with Sand = 30%
(A18b)
Appendix B
The complete site soil survey data can be accessed from the National Cooperative Soil Survey (NCSS; https://ncsslabdatamart.sc.egov.usda.gov/). To access a report, follow these steps:
At the above URL, Enter the “Lab Pedon Number” from table B1 and click “Execute Query” ;
Select the “Lab Pedon Number” 8XP0YYY and Generate Report. Click “Continue” ;
Select “Primary Characterization Report” or other report desired and click “Get Report”.
Table B1. Site number, soil series, location, texture, and National Cooperative Soil Survey Lab Pedon Number (https://ncsslabdatamart.?sc.egov.usda.gov/).
No. Soil series Location Texture Pedon Number
1 Clarion Ames, IA SaL 87P0544
2 Monona Castana, IA SiL 87P0540
3 Cecil 1^[a] Watkinsville, GA SaL 88P0049
4 Sharpsburg Lincoln, NE SiCL 87P0103
5 Hersh Ord, NE SaL 87P0441
6 Keith Albin, WY SiL 87P0263
7 Amarillo Big Spring, TX SaL 87P0455
8 Woodward Woodward, OK L 87P0435
9 Heiden Waco, TX C 87P0647
10 Whitney Fresno, CA SaL 87P0489
11 Academy Fresno, CA SaL 87P0630
12 Los Banos Los Banos, CA SiL 87P0492
13 Portneuf Twin Falls, ID SiL 87P0597
14 Nansene Dusty, WA SiL 87P0497
15 Palouse Pullman, WA SiL 87P0496
16 Zahl Bainville, MT L 87P0532
17 Pierre Wall, SD SiC 87P0525
18 Williams McClusky, ND L 87P0561
19 Barnes Goodrich, ND L 87P0566
20 Sverdrup Wall Lake, MN SaL 87P0576
21 Barnes Morris, MN L 87P0571
22 Mexico Columbia, MO SiL 88P0646
23 Grenada Como, MI SiL 88P0251
24 Tifton Tifton, GA Sa 88P0519
25 Bonifay Tifton, GA Sa 88P0522
26 Cecil 2^[a] Watkinsville, GA SaL 88P0508
27 Hiwassee Watkinsville, GA SaL 88P0515
28 Gaston Salisbury, NC CL 88P0524 & 528
29 Opequon Flintstone, MD CL 88P0503
30 Frederick Hancock, MD C 88P0498
31 Manor Ellicott City, MD L 88P0493
32 Caribou Presque Isle, ME L 88P0722
33 Collamer Ithaca, NY SiL 88P0713
34 Miamian Dayton, OH L 88P0676
35 Lewisburg Columbia, IN CL 88P0147
36 Miami Waveland, IN SiL 88P0152
^[a] Cecil1 was “uneroded” and Cecil2 was “eroded”
To obtain the soil properties of the four satellite pits evaluated on each site, enter the Soil Series name on the NCSS Welcome page, and then select from the list of surveys with similar Lab Pedon Number formats to those in
table B1.
References
Alberts, E. E., Nearing, M. A., Weltz, M. A., Risse, L. M., Pierson, F. B., Zhang, X. C., . . . Simanton, J. R. (1995). Soil component, Chapter 7. In D. C. Flanagan, & M. A. Nearing (Eds.), USDA Water Erosion Prediction Project: Hillslope profile and watershed model documentation. NSERL Report No. 10. West Lafayette, IN: USDA-ARS National Soil Erosion Research Laboratory. Retrieved from https://www.ars.usda.gov/ARSUserFiles/50201000/WEPP/chap7.pdf
Al-Hamdan, O. Z., Pierson, F. B., Nearing, M. A., Williams, C. J., Hernandez, M., Boll, J., . . . Spaeth, K. (2017). Developing a parameterization approach for soil erodibility for the Rangeland Hydrology and Erosion Model (RHEM). Trans. ASABE, 60(1), 85-94. doi:https://doi.org/10.13031/trans.11559
Al-Hamdan, O. Z., Pierson, F. B., Robichaud, P., Elliot, W. J., & Williams, C. J. (2022). New erodibility parameterization for applying WEPP on rangelands using ERMiT. J. ASABE, 65(2), 251-264. doi:https://doi.org/10.13031/ja.14564
Bajracharya, R. M., Elliot, W. J., & Lal, R. (1992). Interrill erodibility of some Ohio soils based on field rainfall simulation. Soil Sci. Soc. Am. J., 56(1), 267-272. doi:https://doi.org/10.2136/sssaj1992.03615995005600010041x
Barnett, A. P., & Rogers, J. S. (1966). Soil physical properties related to runoff and erosion from artificial rainfall. Trans. ASAE, 9(1), 123-128. doi:https://doi.org/10.13031/2013.39895
Bennett H., H. (1939). Soil conservation. New York, NY: McGraw-Hill.
Bjorneberg, D. L. (2001). Evaluating WEPP-predicted furrow irrigation erosion. In J. C. Ascough, & D. C. Flanagan (Ed.), Proc. Soil Erosion Research for the 21st Century (pp. 673-676). St. Joseph, MI: ASABE. doi:https://doi.org/10.13031/2013.4867
Bradford, J. M., Truman, C. C., & Huang, C. (1992). Comparison of three measures of resistance of soil surface seals to raindrop splash. Soil Technol., 5(1), 47-56. doi:https://doi.org/10.1016/0933-3630(92)90006-M
Brazier, R. E., Beven, K. J., Freer, J., & Rowan, J. S. (2000). Equifinality and uncertainty in physically based soil erosion models: Application of the GLUE methodology to WEPP – the Water Erosion Prediction Project – for sites in the UK and USA. Earth Surf. Process. Landforms, 25(8), 825-845. doi:https://doi.org/10.1002/1096-9837(200008)25:8<825::AID-ESP101>3.0.CO;2-3
Buol, S. W., Southard, R. J., Graham, R. C., & McDaniel, P. A. (2003). Soil genesis and classification. Ames, IA: Blackwell.
Dimoyiannis, D. G., Valmis, S., & Vyrlas, P. (2001). A rainfall simulation study of erosion of some calcareous soils. Global Nest: Int. J., 3(3), 179-183.
Douglas-Mankin, K. R., Roy, S. K., Sheshukov, A. Y., Biswas, A., Gharabaghi, B., Binns, A., . . . Daggupati, P. (2020). A comprehensive review of ephemeral gully erosion models. CATENA, 195, 104901. doi:https://doi.org/10.1016/j.catena.2020.104901
Elliot, W. J. (1988). A process-based rill erosion model. PhD Diss. Ames, IA: Iowa State University, Deptartment of Agriculture Engineering. Retrieved from https://dr.lib.iastate.edu/server/api/core/bitstreams/cd5cccac-451e-4587-a32f-2842d3e629b7/content
Elliot, W. J. (2004). WEPP internet interfaces for forest erosion prediction. JAWRA, 40(2), 299-309. doi:https://doi.org/10.1111/j.1752-1688.2004.tb01030.x
Elliot, W. J. (2017). Final report on the Southern Nevada Public Land Management Act Project Number P084: Development of an Online Watershed Interface to predict the effects of forest and fire management on sediment and nutrient loads in surface runoff in the Tahoe Basin. USDA Forest Service, Rocky Mountain Research Station. Retrieved from https://www.fs.usda.gov/psw/partnerships/tahoescience/documents/p084_FinalReport.pdf
Elliot, W. J., Kohl, K. D., & Laflen, J. M. (1988). Methods of collecting WEPP soil erodibility data, ASAE Paper No. MCR:88-138. Proc. Mid-Central Region Mtg. of the ASAE. ASAE.
Elliot, W. J., Laflen, J. M., & Foster, G. R. (1993). Soil erodibility nomographs for the WEPP model, ASAE Paper No. 932046. Proc. 1993 Intl. Summer Mtg. of the ASAE and CSAE. ASABE.
Elliot, W. J., Laflen, J. M., & Kohl, K. D. (1989a). Effect of soil properties on soil erodibility, Paper No. 892150. Proc. 1980 Intl. Summer Mtg. jointly sponsored by the ASAE and the CSAE. ASAE.
Elliot, W. J., Laflen, J. M., Thomas, A. W., & Kohl, K. D. (1997). Photogrammetric and rillmeter techniques for hydraulic measurement in soil erosion studies. Trans. ASAE, 40(1), 157-165. doi:https://doi.org/10.13031/2013.21261
Elliot, W. J., Liebenow, A. M., Laflen, J. M., & Kohl, K. D. (1989b). Compendium of soil erodibility data from WEPP cropland field erodibility experiments 1987 & 1988: NSERL Report No. 3. West Lafayette, IN: USDA-ARS National Soil Erosion Research Laboratory. Retrieved from https://milford.nserl.purdue.edu/weppdocs/comperod/
Elliot, W. J., Olivieri, L. J., Laflen, J. M., & Kohl, K. D. (1990a). Predicting soil erodibility from soil strength measurements. Paper No. 902009. Proc. 1990 Intl. Summer Mtg. of the ASAE. ASAE.
Elliot, W. J., Olivieri, L. J., Laflen, J. M., & Kohl, K. D. (1990b). Predicting soil erodibility from soil properties including classification, mineralogy, climate, and topography. Paper No. 902557. Proc. 1990 Intl. Winter Mtg. of the ASAE. ASAE.
Ellison, W. D. (1947). Soil erosion studies, part I. Agric. Eng., 28(4), 145-146.
El-Swaify, S. A., Dangler, E. W., & Armstrong, C. L. (1982). Soil erosion by water in the tropics. Honolulu, HI: College of Tropical Agriculture and Human Resources, University of Hawaii.
Eriksson, K. A., & Daniels, W. L. (2021). Environmental implications of regional geology and coal mining in the Appalachians. In C. E. Zipper, & J. Skousen (Eds.), Appalachia's coal-mined landscapes: Resources and communities in a new energy era (pp. 27-53). Cham, Switzerland: Springer. doi:https://doi.org/10.1007/978-3-030-57780-3_2
Flanagan, D. C., & Livingston, S. J. (1995). WEPP User Summary Water Erosion Prediction Project, NSERL Report No. 11. West Lafayette, IN: USDA-ARS National Soil Erosion Research Laboratory.
Flanagan, D. C., & Nearing, M. A. (1995). Water Erosion Prediction Project hillslope profile and watershed model documentation. NSERL Report No. 10. West Lafayette, IN: USDA-ARS National Soil Erosion Research Laboratory.
Flanagan, D. C., & Nearing, M. A. (2000). Sediment particle sorting on hillslope profiles in the WEPP model. Trans. ASAE, 43(3), 573-583. doi:https://doi.org/10.13031/2013.2737
Foltz, R. B., Rhee, H., & Elliot, W. J. (2008). Modeling changes in rill erodibility and critical shear stress on native surface roads. Hydrol. Process., 22(24), 4783-4788. doi:https://doi.org/10.1002/hyp.7092
Foster, G. R., & Lane, L. J. (1987). User requirements USDA-Water Erosion Prediction Project (WEPP) Draft 6.3. NSERL Report No. 1. West Lafayette, IN: USDA-ARS National Soil Erosion Research Laboratory. Retrieved from https://www.ars.usda.gov/ARSUserFiles/50201000/WEPP/docs/wepp-user-requirements.pdf
Foster, G. R., & Meyer, L. D. (1972). A closed-form soil erosion equation for upland areas. In J. E. Shen (Ed.), Sedimentation: Symposium to honor Prof. H. A. Einstein. Fort Collins, CO: Colorado State University.
Foster, G. R., Flanagan, D. C., Nearing, M. A., Lane, L. J., Risse, L. M., & Finkner, S. C. (1995). Hillslope erosion component - Chapter 11. In D. C. Flanagan, & M. A. Nearing (Eds.), USDA Water Erosion Prediction Project: Hillslope profile and watershed model documentation. NSERL Report No. 10. West Lafayette, IN: USDA-ARS National Soil Erosion Research Laboratory.
Gilley, J. E., Elliot, W. J., Laflen, J. M., & Simanton, J. R. (1993). Critical shear stress and critical flow rates for initiation of rilling. J. Hydrol., 142(1), 251-271. doi:https://doi.org/10.1016/0022-1694(93)90013-Y
Gui, Y., Zhang, Q., Qin, X., & Wang, J. (2021). Influence of organic matter content on engineering properties of clays. Adv. Civil Eng., 2021, 6654121. doi:https://doi.org/10.1155/2021/6654121
Holz, D. J., Williard, K. W., Edwards, P. J., & Schoonover, J. E. (2015). Soil erosion in humid regions: A review. J. Contemp. Water Res. Educ., 154(1), 48-59. doi:https://doi.org/10.1111/j.1936-704X.2015.03187.x
Huffman, R. L., Fangmeier, D. D., Elliot, W. J., & Workman, S. R. (2013). Chapter 7: Soil erosion by water. In Soil and water conservation engineering (7th ed., pp. 145-170). St. Joseph, MI: ASABE. doi:https://doi.org/10.13031/swce.2013.7
Kinnell, P. I. (1991). The effect of flow depth on sediment transport induced by raindrops impacting shallow flows. Trans. ASAE, 34(1), 161-168. doi:https://doi.org/10.13031/2013.31639
Knisel, W. G. (1980). CREAMS: A field-scale model for Chemicals, Runoff, and Erosion from Agricultural Management Systems, Conservation Research Report No 26. Washington, D.C.: Department of Agriculture, Science and Education Administration.
Laflen, J. M., Elliot, W. J., Simanton, J. R., Holzhey, C. S., & Kohl, K. D. (1991). WEPP: Soil erodibility experiments for rangeland and cropland soils. J. Soil Water Conserv., 46(1), 39-44. Retrieved from https://www.jswconline.org/content/jswc/46/1/39.full.pdf
Lal, R., & Elliot, W. (1994). Soil erosion research methods. In R. Lal (Ed.), Erodibility and erosivity: Chapter 8 (2nd ed.). Ankeny, IA: Soil and Water Conservation Society. doi:https://doi.org/10.1201/9780203739358
Lambe, T. W., & Whitman, R. V. (1969). Soil Mechanics. New York: John Wiley & Sons.
Latifi, N., Rashid, A. S., Siddiqua, S., & Horpibulsuk, S. (2015). Micro-structural analysis of strength development in low- and high swelling clays stabilized with magnesium chloride solution — A green soil stabilizer. Appl. Clay Sci., 118, 195-206. doi:https://doi.org/10.1016/j.clay.2015.10.001
Le Bissonnais, Y. (1996). Aggregate stability and assessment of soil crustability and erodibility: I. Theory and methodology. Eur. J. Soil Sci., 47(4), 425-437. doi:https://doi.org/10.1111/j.1365-2389.1996.tb01843.x
Liebenow, A. M., Elliot, W. J., Laflen, J. M., & Kohl, K. D. (1990). Interrill erodibility: Collection and analysis of data from cropland soils. Trans. ASAE, 33(6), 1882-1888. doi:https://doi.org/10.13031/2013.31553
Liu, X. (2012). Effects of surface coal mining on soil hydraulic properties, Rosebud Mine, Eastern Montana. MS Thesis. Pullman, WA: Washington State University. Retrieved from https://rex.libraries.wsu.edu/esploro/outputs/graduate/Effects-of-surface-coal-mining-on/99900525092201842#file-0
Ma, G., Li, G., Mu, X., Hou, W., Ren, Y., & Yang, M. (2022). Effect of raindrop splashes on topsoil structure and infiltration characteristics. CATENA, 212, 106040. doi:https://doi.org/10.1016/j.catena.2022.106040
Meyer, L. D., & Harmon, W. C. (1979). Multiple-intensity rainfall simulator for erosion research on row sideslopes. Trans. ASAE, 22(1), 100-103. doi:https://doi.org/10.13031/2013.34973
Meyer, L. D., & Harmon, W. C. (1984). Susceptibility of agricultural soils to interrill erosion. Soil Sci. Soc. Am. J., 48(5), 1152-1157. doi:https://doi.org/10.2136/sssaj1984.03615995004800050040x
Meyer, L. D., & McCune, D. L. (1958). Rainfall simulator for runoff plots. Agric. Eng., 39(10), 644-648.
Meyer, L. D., Line, D. E., & Harmon, W. C. (1992). Size characteristics of sediment from agricultural soils. J. Soil Water Conserv., 47(1), 107-111.
Mirzaee, S., & Ghorbani-Dashtaki, S. (2021). Calibrating the WEPP model to predict soil loss for some calcareous soils. Arabian J. Geosci., 14(21), 2198. doi:https://doi.org/10.1007/s12517-021-08646-3
Mohammed, D., & Kohl, R. A. (1987). Infiltration response to kinetic energy. Trans. ASAE, 30(1), 108-111. doi:https://doi.org/10.13031/2013.30410
Nearing, M. A., Foster, G. R., Lane, L. J., & Finkner, S. C. (1989). A process-based soil erosion model for USDA-Water Erosion Prediction Project technology. Trans. ASAE, 32(5), 1587-1593. doi:https://doi.org/10.13031/2013.31195
Nearing, M. A., Wei, H., Stone, J. J., Pierson, F. B., Spaeth, K. E., Weltz, M. A., . . . Hernandez, M. (2011). A rangeland hydrology and erosion model. Trans. ASABE, 54(3), 901-908. doi:https://doi.org/10.13031/2013.37115
Onstad, C. A., & Young, R. A. (1982). Erosion characteristics of three northwestern soils. Trans. ASAE, 25(2), 367-371. doi:https://doi.org/10.13031/2013.33537
Pierson, F. B., Robichaud, P. R., Moffet, C. A., Spaeth, K. E., Hardegree, S. P., Clark, P. E., & Williams, C. J. (2008). Fire effects on rangeland hydrology and erosion in a steep sagebrush-dominated landscape. Hydrol. Process., 22(16), 2916-2929. doi:https://doi.org/10.1002/hyp.6904
Rahimi, H., Pazira, E., & Tajik, F. (2000). Effect of soil organic matter, electrical conductivity and sodium adsorption ratio on tensile strength of aggregates. Soil Tillage Res., 54(3), 145-153. doi:https://doi.org/10.1016/S0167-1987(00)00086-6
Rawls, W. J., Brakensiek, D. L., Simanton, J. R., & Kohl, K. D. (1990). Development of a crust factor for a Green Ampt model. Trans. ASAE, 33(4), 1224-1228. doi:https://doi.org/10.13031/2013.31461
Renard, K. G., Foster, G. R., Weesies, G. A., McCool, D. K., & Yoder, D. C. (1997). Predicting soil erosion by water: A guide to conservation planning with the Revised Universal Soil Loss Equation (RUSLE). Agricultural Handbook No. 703. Washington, D.C.: USDA. Retrieved from https://www.ars.usda.gov/ARSUserFiles/64080530/RUSLE/AH_703.pdf
Robichaud, P. R., Elliot, W. J., Pierson, F. B., Hall, D. E., & Moffet, C. A. (2007). Predicting postfire erosion and mitigation effectiveness with a web-based probabilistic erosion model. CATENA, 71(2), 229-241. doi:https://doi.org/10.1016/j.catena.2007.03.003
Robichaud, P. R., Wagenbrenner, J. W., & Brown, R. E. (2010). Rill erosion in natural and disturbed forests: 1. Measurements. Water Resour. Res., 46(10). doi:https://doi.org/10.1029/2009WR008314
SAS Institute Inc. (1982). SAS user's guide, Statistics, 1982 Edition. Cary, NC: SAS Institute.
Savabi, M. R., & Williams, J. R. (1995). Water balance and percolation - Chapter 5. In D. C. Flanagan, & M. A. Nearing (Eds.), USDA Water Erosion Prediction Project: Hillslope profile and watershed model documentation. NSERL Report No. 10. West Lafayette, IN: USDA-ARS National Soil Erosion Research Laboratory.
Smerdon, E. T., & Beasley, R. P. (1961). Critical tractive forces in cohesive soils. Agric. Eng., 42(1), 26-29.
Snedecor, G. W., & Cochran, W. G. (1972). Statistical Methods (6th ed.). Ames, IA: Iowa State University.
Swanson, N. P. (1965). Rotating-boom rainfall simulator. Trans. ASAE, 8(1), 71-72. doi:https://doi.org/10.13031/2013.40430
Tiscareno-Lopez, M., Weltz, M. A., & Lopes, V. L. (1995). Assessing uncertainties in WEPP's soil erosion predictions on rangelands. J. Soil Water Conserv., 50(5), 512-516. Retrieved from https://www.jswconline.org/content/jswc/50/5/512.full.pdf
USDA Agricultural Research Service (USDA-ARS). (2008). Draft user's reference guide, Revised Universal Soil Loss Equation version 2 (RUSLE2). Retrieved from https://www.nrcs.usda.gov/resources/education-and-teaching-materials/soil-facts
USDA Natural Resource Conservation Service (NRCS). (2022a). Official soil series descriptions (OSD) fact sheet. Retrieved from https://www.nrcs.usda.gov/wps/portal/nrcs/detail/soils/survey/class/data/?cid=nrcs142p2_053586
USDA Natural Resource Conservation Service (NRCS). (2022b). Soil formation and classification. Retrieved from https://www.nrcs.usda.gov/wps/portal/nrcs/detail/soils/edu/?cid=nrcs142p2_054278
USDA Soil Conservation Service (USDA-SCS). (1990). Water Erosion Prediction Project soil characterization data, Volumes 1-8. Lincoln, NE: USDA Soil Conservation Service National Soil Survey Laboratory. Retrieved from https://ncsslabdatamart.sc.egov.usda.gov
Wischmeier, W. H., & Mannering, J. V. (1969). Relation of soil properties to its erodibility. Soil Sci. Soc. Am. J., 33(1), 131-137. doi:https://doi.org/10.2136/sssaj1969.03615995003300010035x
Wischmeier, W. H., & Smith, D. D. (1965). Predicting rainfall erosion losses from cropland east of the Rocky Mountains: Guide for selection of practices for soil and water conservation. Agriculture Handbook No. 282. Washington, D.C.: USDA-Agricultural Research Service.
Wischmeier, W. H., & Smith, D. D. (1978). Predicting rainfall erosion losses: A guide to conservation planning. Agriculture Handbook No 537. Washington, D.C.: U.S. Deptartment of Agriculture. Retrieved from https://naldc.nal.usda.gov/catalog/CAT79706928
Yalin, Y. M. (1963). An expression for bed-load transportation. J. Hydraul. Div., 89(3), 221-250. doi:https://doi.org/10.1061/JYCEAJ.0000874
Yao, C., Lei, T., Elliot, W. J., McCool, D. K., Zhao, J., & Chen, S. (2008). Critical conditions for rill initiation. Trans. ASABE, 51(1), 107-114. doi:https://doi.org/10.13031/2013.24231

American Society of
Agricultural and Biological Engineers

2950 Niles Road, St. Joseph, MI 49085
Phone: +12694290300 Fax: +12694293852
Copyright © 2025 American Society of Agricultural and Biological Engineers

Table 1. Site number, soil series, location, erodibility, and textural properties of WEPP cropland soils. Definitions of variables and their units are presented in table 2.
No.	Soil series	Location	Texture	K_i1^[a]	K_i2	K_r	t_c	Clay	Silt	Sand	FSi	VfSa	MSa	Rock	OC
1	Clarion	Ames, IA	FSL	2.03	2.03	4.6	0.4	18.9	27.2	53.9	14.4	10.9	13.8	4	1.05
2	Monona	Castana, IA	SiL	1.78	1.78	7.6	2.8	20.1	74.8	5.1	27	4.9	0	0	1.02
3	Cecil 1^[b]	Watkinsville, GA	SaL	1.18	1.18	8.4	2.2	8.6	16.7	74.7	11.5	5.9	22.4	4	0.55
4	Sharpsburg	Lincoln, NE	SiCL	1.65	1.85	5.3	3.2	39.8	55.4	4.8	23.9	4.6	0	0	1.85
5	Hersh	Ord, NE	SaL	3.83	3.93	11.2	1.7	9.6	13.4	77	2.2	32.9	18.2	0	0.49
6	Keith	Albin, WY	SiL	3.37				19.3	31.8	48.9	11.3	44	1	0	0.91
7	Amarillo	Big Spring, TX	SaL	3.82	4.12	45.3	1.7	7.3	7.7	85	2.9	21.1	23.9	0	0.16
8	Woodward	Woodward, OK	L	3.88	4	25	1.3	12.3	39.9	47.8	9.1	39	1.1	0.3	0.82
9	Heiden	Waco, TX	C	2	1.7	8.9	2.9	53.1	38.3	8.6	29.3	4.5	1	0	1.36
10	Whitney	Fresno, CA	SaL	2.66	2.74	23.3	4.7	7.2	21.7	71.1	9.7	8.1	19.2	4	0.19
11	Academy	Fresno, CA	SaL	3.06	2.88	5.7	1.6	8.2	29.1	62.7	12.2	20.2	13.7	4	0.41
12	Los Banos	Los Banos, CA	SiL	2.53	2.5	0.6	2.1	43	41.0	16	21.1	11.3	0.8	0	1.47
13	Portneuf	Twin Falls, ID	SiL	1.26	1.26	10.6	3.1	11.1	67.4	21.5	29.4	19.3	0.5	0	0.72
14	Nansene	Dusty, WA	SiL	3.13	3.12	30.7	3.1	11.1	68.8	20.1	30.6	18.1	0.3	0	1.49
15	Palouse	Pullman, WA	SiL	3.74	4.32	6.6	0.7	20.1	70.1	9.8	30.9	8.8	0.1	0	1.76
16	Zahl	Bainville, MT	L	3.21	3.17	12.3	3.5	24	29.7	46.3	14.7	12.5	7.2	9	1.69
17	Pierre	Wall, SD	SiC	2.17	2.18	11.7	4.8	49.5	40.9	9.6	21.4	7.3	0.8	0	1.46
18	Williams	McClusky, ND	L	2.95	2.94	4.5	3.4	26	32.4	41.6	17.5	11.5	8.7	5	1.79
19	Barnes	Goodrich, ND	L	1.72	1.71	3.3	2.3	24.6	36.0	39.5	20.6	12.1	6.7	4	3.26
20	Sverdrup	Wall Lake, MN	SaL	2.14	2.11	10	1.4	7.9	16.8	75.3	9.4	3.7	27.3	0.3	1.28
21	Barnes	Morris, MN	L	1.6	1.6	6.3	4	17	34.4	48.6	16.5	11.4	10.2	6	1.98
22	Mexico	Columbia, MO	SiL	2.97	2.97	3.6	0.7	26	68.7	5.3	33.1	1.1	1.5	0	1.56
23	Grenada	Como, MS	SiL	2.64	2.63	7.3	4.5	20.2	77.8	2	36.4	1.5	0.2	0	1.27
24	Tifton	Tifton, GA	Sa	0.76	0.77	11.3	3.5	2.8	10.8	86.4	4.8	13.3	26.2	23	0.46
25	Bonifay	Tifton, GA	Sa	1.58	0.87	17.9	1	3.3	5.5	91.2	3.1	16.2	18.6	1	0.32
26	Cecil 2^[b]	Watkinsville, GA	SaL	1.92	1.86	3.8	4.5	19.8	15.6	64.6	9.3	5.9	18.5	6	0.70
27	Hiwassee	Watkinsville, GA	SaL	1.88	1.88	10.3	2.3	14.7	21.6	63.7	15.3	4.3	19.2	3	0.83
28	Gaston	Salisbury, NC	CL	2.19	2.04	4.9	4.4	39.1	25.4	35.5	17.6	7.5	10.1	0.3	1.12
29	Opequon	Flintstone, MD	CL	3.3	3.2	3.5	6.3	31.1	31.2	37.7	22.8	5.9	9.6	14	1.42
30	Frederick	Hancock, MD	C	2.43	2.48	8.4	6.6	16.6	58.3	25.1	39.4	5.2	6.6	14	1.32
31	Manor	Ellicot City, MD	L	2.68	2.69	5.4	3.6	25.7	30.7	43.6	23.4	7.1	10.6	8	0.96
32	Caribou	Presque Isle, ME	L	1.55	1.45	4.5	4.3	12.2	40.8	47	25.6	11.5	7.8	47	2.28
33	Collamer	Ithaca, NY	SiL	3.74	3.46	24.13	6.4	15	78.0	7	47.5	4.6	0.7	0.3	1.01
34	Miamian	Dayton, OH	L	1.59	1.65	9.6	5.5	25.3	44.1	30.6	28.6	6.4	7.4	3	1.75
35	Lewisburg	Columbia, IN	CL	2.42	2.47	5.9	3.4	29.3	32.2	38.5	19.4	10.9	7.2	6	0.87
36	Miami	Waveland, IN	SiL	2.03	1.97	9.5	3.3	23.1	72.7	4.2	36.9	2	0.8	0	0.82
^[a] Units for all properties are presented in table 2 ^[b]“Cecil 1” was a “non-eroded” site, and “Cecil 2” was an “eroded” site, but both soils were classified as “Cecil”

Table 2. Definitions of soil properties and their units for the 1986-88 WEPP cropland erodibility study that were selected by at least one of the stepwise procedures for inclusion in regression equations to estimate soil erodibility in Elliot et al. (1988, 1989a, 1990a, b, 1993).
Variable	Definition	Units
15Bar	Water content of soil at 15 Bar suction pressure often considered a measure of wilting point soil water content for the fine fraction (<2 mm dia)	Percent
AgStab	Aggregate Stability for aggregates < 5 mm	Percent
Al	Aluminum content in fine fraction	Percent
Ca	Calcium content in fine fraction	MEQ/100 g
CaCO₃	Calcium carbonate content in fine fraction	Percent
CEC	Cation exchange capacity	MEQ/100 g
Clay	Clay content (<0.002 mm dia) of fine fraction (<2 mm dia)	Percent
Cond	Electrical conductivity of fine fraction	mmhos/cm
Depth	Depth of A Horizon	cm
Fe	Iron content in fine fraction	Percent
FSi	Fine silt content (0.002 mm – 0.02 mm dia) of fine fraction (<2 mm dia)	Percent
K or K-Factor	Soil erodibility value used in the USLE and RUSLE models. English units used in analyses (ton acre h [hundreds of acre-ft tonf in.]^-1)
K_i1	Interrill erodibility prior to 1989 x 10^-6	kg s m^-4
K_i2	Interrill erodibility from 1989 to present x 10^-6	kg s m^-4
K_r	Rill erodibility x 10³	s m^-1
LL	Liquid limit water content for soil fraction < 0.4 mm	Percent
M	Wischmeier and Smith’s (1978) texture property = (Silt + VfSa) x (100-Clay)	Percent²
MC	Soil water content measured in the field prior to the start of rainfall simulation	g water/g soil
Mg	Magnesium content in fine fraction	MEQ Mg/100 g
Miner	Mineralogy of the dominant clay fraction: 1 – Kaolinitic; 2 – Mixed; 3 – Mixed-Smectitic; 4 – Smectitic and Montmorillinitic; 5 – Calcitic (Elliot et al., 1990b)
MSa	Medium sand content (0.25 mm - 0.5 mm dia) of fine fraction (<2 mm dia)	Percent
N	Nitrogen content in fine fraction	Percent
Na	Sodium content in fine fraction	MEQ Na/100 g
OC	Organic Carbon in fine fraction (< 2 mm dia)	Percent
OM	Organic matter in fine fraction (< 2 mm dia)	Percent
Order	Taxonomic Order of soil: 1 – Vertisols and Aridisols; 2 – Ultisols; 3 – Mollisols, Inceptisols, Alfisols, Entisols and Spodosols (Elliot et al., 1990b)
PL	Plastic limit water content for soil fraction < 0.4 mm	Percent
PI	Plasticity Index = LL - PL	Percent
Rock	Rock content of whole soil sample (> 2 mm dia)	Percent
Sand	Sand content (0.05 mm – 2 mm dia) of fine fraction (<2 mm dia)	Percent
SAR	Sodium Adsorption Ratio (Na/((Ca+Mg)/2)^0.5); Measured on soils 1, 2, 4-22, and 33. Otherwise set to default value of 0.1	-
Silt	Silt content (0.002 mm – 0.05 mm dia) of fine fraction (<2 mm dia)	Percent
Slope	Steepness of the study site.	Percent
SpSf	Specific surface of soil particles measured with ethylene glycol monoethyl ether (EGME)	mg g^-1 EGME
t_c	Critical shear of soil	Pa
VSJ	Handheld vane shear measurement of soil strength on field plot	Pa
VfSa	Very fine sand content (0.05 mm – 0.01 mm dia) of fine fraction (<2 mm dia)	Percent
WDClay	Water dispersible clay in fine fraction	Percent
WDSilt	Water dispersible silt in fine fraction	Percent

Table 3a. Physical properties that were useful for predicting WEPP cropland soil erodibility values. The definitions of the variables and their units are presented in table 2. The taxonomic orders that were split for estimating K_i2 in equation 4a are highlighted.
No	M	WD Clay	WDSilt	PL	LL	PI	15Bar	Ag Stab	Miner	Miner/Clay	MC	Slope	Order	VSJ
1	3090	3.1	38.1	18	29	11	8.8	7	3	0.159		5.17	Hapludoll
2	6368	2.9	58.4	23	37	14	10.2	9	3	0.149		12.62	Hapludoll
3	2066	7.5	19.3	15	18	3	3.6	10	1	0.116		4	Kanhapludult
4	3612	7.5	62.7	25	19	11	17.4	4	4	0.101	0.06	5.68	Argiudoll	51972
5	4186	1.1	18.4	17	37	14	4.8	9	3	0.313	0.067	6.61	Ustorthent	7844
6	6117	3.4	40.7	19	18	3	10.3	6	3	0.155	0.03	5.08	Argiustoll	27456
7	2670	3.5	10.6	17	40	15	3.4	3	2	0.274	0.11	3.59	Paleustalf	13238
8	6920	3.2	43	24	18	1	5.9	3	3	0.244	0.247	7.14	Haplustept	21573
9	2007	17.6	59.6	20	29	10	18.9	61	4	0.075	0.003	3.86	Chromustert	19121
10	2765	4.1	23.4	13	18	1	2.6	13	2	0.278	0	7.36	Haploxeralf	52167
11	4526	2.6	34.2	17	28	4	3.5	9	2	0.244	0.033	4.51	Haploxeralf	31379
12	2981	12.3	54.5	21	55	35	19.1	15	4	0.093	0.03	4.03	Haploxeralf	68543
13	7708	3.5	77.6	23	15	2	8.3	33	2	0.180	0.053	5.57	Haplocalcid	67661
14	7725	2.2	71.7	24	18	1	7.3	20	2	0.180	0.097	6.1	Haploxeroll	65209
15	6304	2	69.8	23	49	28	9	20	2	0.100	0.13	6.47	Haploxeroll	38243
16	3207	6.2	36.9	18	27	4	10.2	12	3	0.125	0.133	7.56	Calciustoll	56875
17	2434	14.3	53.3	14	28	4	19.4	40	4	0.081	0.18	6.65	Haplustept	45304
18	3249	8.8	44.8	19	32	9	11.9	5	3	0.115	0.15	5.05	Argiustoll	73545
19	3627	7.1	43.8	27	34	16	12.7	13	3	0.122	0.09	5.78	Hapludoll	55894
20	1888	2.8	21.7	23	24	10	5.4	13	2	0.253	0.093	4.23	Hapludoll	57855
21	3801	8.6	41.5	18	36	17	9.3	6	3	0.176	0.17	8.27	Hapludoll	55894
22	5165	6.4	82.8	23	40	13	11.5	6	3	0.115	0.1	3.87	Epiaqualf	17651
23	6328	5	89.9	23	33	10	10.9	5	3	0.149	0.01	8.68	Fraglossudalf	17651
24	2343	2.5	11.8	17	34	16	2.1	3	1	0.357,	0.065	4.53	Kandiudult	33340
25	2098	1.5	7.6	17	40	17	1.1	3	1	0.303	0.03	3.9	Paleudult	37262
26	1724	14.4	18.4	17	33	10	8.3	15	1	0.051	0.01	4.32	Kanhapludult	56875
27	2209	10.6	26.4	16	18	1	6.4	3	1	0.068	0.11	3.96	Kanhapludult	33340
28	2004	24.5	32	20	18	1	14.8	23	1	0.026	0.115	6.38	Hapludult	39224
29	2556	24.4	58	26	29	12	14.2	29	2	0.064	0.13	12.05	Hapludalf	29418
30	5296	8.4	60	26	25	9	6.8	18	2	0.120	0.115	12.85	Paleudult	37263
31	2809	12.6	41.6	23	39	19	11.7	6	2	0.078	0.065	8.54	Dystrudept	3137
32	4592	7.7	41.9	29	38	12	8.9	10	2	0.164	0.135	8.85	Haplorthod	23534
33	7021	6.9	81.7	23	30	8	6.9	10	2	0.133	0.025	8.68	Hapludalf	23534
34	3772	12.2	52.9	21	44	21	10.9	10	2	0.079	0.11	8.92	Hapludalf	31379
35	3047	19.4	40.6	16	32	3	11.3	18	1	0.034	0.025	7.45	Hapludalf	15689
36	5744	7.5	86.6	21	29	6	10	5	3	0.130		6.5	Hapludalf	23534

Table 3b. Chemical properties for soils presented in table 1 that were useful for predicting soil erodibility in the WEPP cropland soil erodibility study. Definitions and units for the variables are presented in table 2.
No	Fe	Al	N	Ca	Na	Mg	CEC	SAR	Cond	CaCO₃	Sp Sf
1	0.9	0.1	0.106	12.1	0.1	3.6	14.1	0.3	0.48	0.0	19
2	1.2	0.1	0.112	17.2	0.1	4.8	18.9	0.1	0.46	0.0	19
3	0.9	0.1	0.049	2.3	0.03	0.7	3	0.1	0.6	0.0	4
4	1.1	0.2	0.171	19.4	0.1	6.6	29.4	0.3	0.2	0.5	36
5	0.2	0.1	0.057	5.6	0	1.5	7.6	0.3	0.04	0.5	7
6	0.5	0.1	0.093	12.9	0.03	3.1	18.3	0.1	0.38	0.5	15
7	0.2	0.03	0.045	5.9	0.3	0.5	5.1	0.3	0.25	0.5	6
8	0.6	0.03	0.092	0	0.03	1.9	10.4	0.3	0.56	1.0	8
9	0.6	0.1	0.125	0	0.1	1	33.3	0.3	0.66	45.	38
10	0.5	0.03	0.024	2.8	0.1	0.7	3.4	1	0.11	0.5	3
11	0.9	0.1	0.034	3.2	0	1.4	5.3	0.3	0.84	0.5	4
12	0.8	0.1	0.146	28.3	0.5	8.6	38.7	1	0.89	1.0	45
13	0.3	0.03	0.083	0	0.2	6.3	12.6	1	0.76	23.	11
14	0.9	0.1	0.112	8.6	0.1	2.3	15.1	0.3	0.93	0.5	11
15	0.9	0.1	0.149	12.7	0.1	2.5	19.6	0.3	1.37	0.5	18
16	0.9	0.1	0.173	22.6	0.1	4.3	20	0.3	0.72	0.5	23
17	1.5	0.1	0.16	32.5	0.1	7.1	36.8	0.3	0.5	0.5	46
18	0.7	0.1	0.173	24.6	0.1	6.2	22.7	0.3	0.81	0.5	25
19	0.9	0.1	0.298	13.8	0.03	5	24.1	0.3	0.42	0.5	21
20	0.5	0.1	0.108	9.5	0	2.1	11	0.3	0.35	0.5	8
21	0.8	0.1	0.181	0	0	5.1	19.5	0.3	0.72	1.0	15
22	1.6	0.2	0.152	15.6	0.1	3	21.3	0.3	0.5	0.5	34
23	1.7	0.2	0.129	4.1	0.1	1.5	12	0.1	0.3	0.5	17
24	0.5	0.1	0.034	1.5	0.03	0.3	2.1	0.1	0.2	0.5	3
25	0.2	0.1	0.015	0.7	0.1	0.1	1.7	0.1	0.2	0.5	3
26	2.4	0.2	0.057	2.2	0.03	0.6	3.5	0.1	0.2	0.5	8
27	1.5	0.2	0.075	1.6	0.03	0.4	4.4	0.1	0.2	0.5	7
28	4.5	0.5	0.107	4.8	0.03	1.8	9.2	0.1	0.2	0.5	18
29	2.7	0.3	0.16	10.3	0.1	1	12.9	0.1	0.2	0.5	18
30	1.4	0.2	0.135	5	0	1.1	8.4	0.1	0.2	0.5	10
31	3.1	0.4	0.108	7.9	0.03	1.7	11.9	0.1	0.2	0.5	17
32	1.3	0.4	0.193	5.5	0.1	1.5	12	0.1	0.2	0.5	15
33	1.1	0.2	0.105	7.8	0.1	0.5	8.9	0.3	0.2	0.5	10
34	1.4	0.1	0.153	10	0	4.5	12.5	0.1	0.2	0.5	21
35	1.6	0.2	0.099	11.2	0.03	3.2	14.4	0.1	0.2	9.0	17
36	1.4	0.2	0.091	11	0.1	1.2	13.3	0.1	0.2	0.5	20

Table 4. WEPP cropland erodibility study soil properties evaluated but not selected in stepwise multiple regression procedures. Correlation coefficients (r) greater than 0.25 are noted.
Elliot et al. (1988)
	Temperature on day of study (r = -0.31 with K_i2)
	Final soil water content; Depth of A Horizon
	Sodium content and density of rainfall and runoff simulation water
Elliot et al. (1989a)
	Fine sand (r = 0.27 with K_r and -0.27 with t_c); Coarse sand (r = -0.31 with K_i2)
	Very coarse sand (r = -0.32 with K_i2); USLE K-Factor (r = 0.357 with K_i2) Soil water content at 1/3 bar (field capacity) (r = -0.32 with K_r and 0.27 with t_c)
Elliot et al. (1990a)
	Coefficient of linear expansion (r = -0.30 with K_r)
	Onsite torvane on rill bottom before rainfall (r = -0.40 with K_r) Onsite torvane on rill side before rainfall (r = -0.41 with K_r) Onsite torvane on interrill plot before rainfall (r = -0.29 with K_r) Onsite torvane on rill bottom after rainfall (r = -0.34 with t_c) Onsite torvane, pocket penetrometer and fall cone penetrometer in saturated soil outside of plots Onsite pocket penetrometer reading of soil strength with big foot before rainfall (r = -0.26 with K_i2)
	Onsite pocket penetrometer reading of soil strength with big foot after rainfall (r = -0.47 with K_i2) Onsite fall cone penetrometer measurement of soil strength (r = -0.52 with K_i2)
	Laboratory unconfined consolidation stress (r = 0.29 with t_c)
	Laboratory soil stress at 1% consolidation (r = -0.36 with K_r) and 4% consolidation (r = -0.33 with K_r and 0.34 with t_c)
	Laboratory pinhole critical shear test (r = -0.36 with K_r); Laboratory angle of internal shear
	Coefficient of Linear Expansion (r = -0.30 with K_r) PI/Clay (r = -0.27 with K_i2, 0.36 with K_r and -0.31 with t_c); Sp Sf / Clay (r = -0.25 with t_c) Clay – WDClay (r = -0.32 with K_r); (Clay – WDClay)/Clay (r = 0.42 with K_i2 and -0.42 with t_c)
Elliot et al. (1990b)
	Average annual temperature of site; Average annual precipitation of site (r = -0.39 with K_i2)

Table 5. Correlation coefficients (r) between K_i1, K_i2, K_r, t_c, and clay content and soil properties in Elliot et al.'s (1989a, 1990a, 1990b, 1993) analyses for WEPP cropland soil erodibility studies. Soil property descriptions and units for variables are described in table 2. The soil properties with the greatest correlation coefficient for each erodibility value are highlighted.
Property	K_i1	K_i2	K_r	t_c	Clay
K_i1	1.000	0.979	0.364	-0.171	-0.031
K_i2	0.979	1.000	0.405	-0.102	-0.034
K_r	0.364	0.405	1.000	0.049	-0.417
t_c	-0.171	-0.102	0.049	1.000	0.196
Clay	-0.031	-0.034	-0.417	0.196	1.000
Silt	0.125	0.166	-0.076	0.209	0.222
Sand	-0.085	-0.118	0.252	-0.256	-0.634
FSi	0.034	0.063	-0.155	0.442	0.350
VfSa	0.415	0.368	0.249	-0.558	-0.327
MSa	-0.210	-0.188	0.201	-0.070	-0.549
M	0.338	0.336	0.168	-0.083	-0.227
USLE-K Factor	0.337	0.357	0.183	-0.008	-0.166
Rock	-0.284	-0.265	-0.189	0.314	-0.198
OC	-0.121	-0.072	-0.408	0.156	0.437
OM	-0.008	0.027	-0.395	0.186	0.410
WD Clay	-0.117	-0.126	-0.350	0.495	0.698
WDSilt	0.151	0.174	-0.162	0.233	0.385
WDClay/Clay	-0.406	-0.418	-0.002	0.451	-0.183
CEC/Clay	0.087	0.123	0.091	-0.360	-0.161
Fe	-0.046	-0.032	-0.355	0.430	0.426
Al	-0.082	-0.099	-0.373	0.373	0.291
N	-0.042	0.009	-0.398	0.194	0.515
Ca	0.160	0.187	-0.256	-0.042	0.563
Na	0.128	0.209	0.196	-0.082	0.244
Mg	-0.149	-0.121	-0.343	-0.061	0.509
CEC	0.019	0.025	-0.338	-0.073	0.808
SAR	0.073	0.127	0.150	-0.052	0.032
Cond	0.184	0.236	-0.075	-0.321	0.119
CaCO₃	-0.180	-0.205	-0.035	0.002	0.365
Ag Stab	-0.070	-0.074	-0.099	0.256	0.579
SpSf	-0.037	-0.028	-0.391	0.045	0.899
15Bar	-0.051	-0.055	-0.477	0.174	0.964
PL	0.062	0.086	-0.171	0.177	0.101
LL	0.002	0.030	-0.174	-0.061	0.147
PI	-0.033	0.026	-0.220	-0.065	0.147
MC	0.160	0.220	0.118	0.059	0.029
Slope	0.103	0.141	-0.096	0.601	0.058
Order	0.486	0.519	0.054	-0.088	0.025
Miner	0.191	0.211	-0.161	-0.151	0.559
Miner/Clay	0.071	0.053	0.490	-0.300	-0.731
VSJ	-0.297	-0.275	-0.116	0.075	0.056

Table 6. Soil properties selected by stepwise procedures and coefficients of determination (r²) for each stepwise multiple regression analysis in WEPP cropland soils erodibility study. Definitions and units for the variables are presented in table 2. Complete equations are in Appendix A.
	Predictive variables to estimate soil erodibility
Analysis	K_i1or K_i2	K_r	t_c
Elliot al., 1988 for soil nos. 4-21	For K_i1: Clay, VSJ; r² = 0.35	Clay, VSJ; r² = 0.35	MC, Depth; r² = 0.16
Elliot et al., 1989a	For K_i1: Clay, AgStab, WDClay, Mg, Fe, Al, Cond; r²=0.73	M, CEC, SAR, Al, OC; r²=0.71	Clay, MC, VfSa, CaCO₃, SAR, SpSf, WDClay, Sand; r²=0.61
Elliot et al., 1990a	For K_i2: CaCO₃, Na, Ca, Mg, PI, Clay, SpSf, WDClay; r²=0.55	M, Sand, LL, Al, Mg; r²=0.76	Cond, Sand, VSJ, WDClay, Clay, Silt; r²=0.59
Elliot et al., 1990b	For K_i2: Miner, Cond, WDSilt, FSilt, CaMg, Na, Clay, WDClay, SpSf; r²=0.79	Miner, Clay, Al, Mg, 15Bar, M, VfSa; r²=0.76	Clay, WDClay, Sloper²=0.75
Elliot et al., 1993	For K_i2: Silt, VfSa, Sand, OC, PL, Mg; r²=0.35	Silt, VfSa, Sand, OC, PL, Mg; r²=0.81	Silt, VfSa, Sand, OC, PL, Mg; r²=0.26
Alberts et al., 1995	For K_i2: Sand, VfSa, Clay; r²=0.24	Sand, VfSa, OM; r²=0.55	Sand, VfSa, Clay r²=0.23

Table B1. Site number, soil series, location, texture, and National Cooperative Soil Survey Lab Pedon Number (https://ncsslabdatamart.?sc.egov.usda.gov/).
No.	Soil series	Location	Texture	Pedon Number
1	Clarion	Ames, IA	SaL	87P0544
2	Monona	Castana, IA	SiL	87P0540
3	Cecil 1^[a]	Watkinsville, GA	SaL	88P0049
4	Sharpsburg	Lincoln, NE	SiCL	87P0103
5	Hersh	Ord, NE	SaL	87P0441
6	Keith	Albin, WY	SiL	87P0263
7	Amarillo	Big Spring, TX	SaL	87P0455
8	Woodward	Woodward, OK	L	87P0435
9	Heiden	Waco, TX	C	87P0647
10	Whitney	Fresno, CA	SaL	87P0489
11	Academy	Fresno, CA	SaL	87P0630
12	Los Banos	Los Banos, CA	SiL	87P0492
13	Portneuf	Twin Falls, ID	SiL	87P0597
14	Nansene	Dusty, WA	SiL	87P0497
15	Palouse	Pullman, WA	SiL	87P0496
16	Zahl	Bainville, MT	L	87P0532
17	Pierre	Wall, SD	SiC	87P0525
18	Williams	McClusky, ND	L	87P0561
19	Barnes	Goodrich, ND	L	87P0566
20	Sverdrup	Wall Lake, MN	SaL	87P0576
21	Barnes	Morris, MN	L	87P0571
22	Mexico	Columbia, MO	SiL	88P0646
23	Grenada	Como, MI	SiL	88P0251
24	Tifton	Tifton, GA	Sa	88P0519
25	Bonifay	Tifton, GA	Sa	88P0522
26	Cecil 2^[a]	Watkinsville, GA	SaL	88P0508
27	Hiwassee	Watkinsville, GA	SaL	88P0515
28	Gaston	Salisbury, NC	CL	88P0524 & 528
29	Opequon	Flintstone, MD	CL	88P0503
30	Frederick	Hancock, MD	C	88P0498
31	Manor	Ellicott City, MD	L	88P0493
32	Caribou	Presque Isle, ME	L	88P0722
33	Collamer	Ithaca, NY	SiL	88P0713
34	Miamian	Dayton, OH	L	88P0676
35	Lewisburg	Columbia, IN	CL	88P0147
36	Miami	Waveland, IN	SiL	88P0152
^[a] Cecil1 was “uneroded” and Cecil2 was “eroded”

Estimating WEPP Cropland Erodibility Values From Soil Properties

Published in Journal of the ASABE 66(2): 329-351 (doi: 10.13031/ja.15218). 2023 American Society of Agricultural and Biological Engineers.

The Water Erosion Prediction Project (WEPP) Model

Soil Erodibility

Methods

Results

Discussion

Conclusions

Acknowledgments

Appendix A

Appendix B

References