Estimating soil erodibility for the RUSLE with rainfall simulation in central Queensland, Australia

Site description

Data were collected from two beef cattle grazing properties within the semi-arid Fitzroy Basin: (1) the Brigalow Catchment Study and surrounds in 2017; and (2) Belmont Research Station in 2021 (Fig. 3).

Fig. 3.

Relative location of the Fitzroy Basin in the Great Barrier Reef catchments, Queensland, Australia (left), and study locations within the Fitzroy Basin (right).

The Brigalow Catchment Study is located near Theodore in central Queensland (24.81°S, 149.80°E, GDA 94) and was described in detail in Cowie et al. (2007). It is located in the Dawson River catchment of the Fitzroy Basin, and the climate is semi-arid to subtropical with a mean annual rainfall of 650 mm. Average maximum monthly temperature in summer is 33°C (December–February) and average minimum temperature in winter of 6.5°C (June–August). Four sites were investigated in this location: two Brown Sodosols and two Black Vertosols. Sites with similar soil types showed enough soil texture variation to be considered separately as Vertosol₁ and Vertosol₂, and Sodosol₁ and Sodosol₂. These soil types are representative of 56% of the Fitzroy Basin with the land use of grazing: Vertosols (28%) and Sodosols (28%) (Roots 2016; Elledge and Thornton 2017). The two Sodosols in this study have a low average exchangeable sodium percentage (ESP) in the 0–0.1 m depth of 0.4% for Sodosol₁ and 4.5% for Sodosol₂. Vertosol₁ showed ESP at the depth of 0–0.1 m was 6.8%, indicating marginal sodicity (Hazelton and Murphy 2007).

Belmont Research Station is located approximately 37 km north of Rockhampton adjacent to the Fitzroy River in the Fitzroy catchment of the Fitzroy Basin (23.23°S, 150.39°E GDA 94) and was described in detail in Coventry and Murtha (1976). The station is located close to the Tropic of Capricorn, and the climate is sub-tropical with a mean annual rainfall of 800 mm. Average maximum monthly temperature in summer is 32°C (December–February), and average minimum monthly temperature in winter is 9°C (June–August). Two sites and soil types were investigated in this location including a Leptic Tenosol and a Black Dermosol. These soils are collectively representative of 14% of the Fitzroy Basin under grazing: Dermosols (11%) and Tenosols (3%) (Roots 2016).

Treatments

All soils were studied in situ to preserve soil consolidation and organic bonds. The plot surface was prepared as an artificially scalded surface to represent consolidated bare ground by mechanically removing the surface organic A1 horizon to expose the hard-set layer beneath and was left to weather and reconsolidate for up to three months. Use of artificial scalds in place of natural scalds was validated in Loch (2015a) and Loch (2017).

The Vertosol and Sodosol soils had an additional natural scald treatment as well as an artificial scald. This treatment identified bare surfaces which had been weathered in situ for at least a year, eroding the remnant A surface horizon to a hard-set layer.

Rainfall simulation, runoff generation, and water quality sampling

Plots were bunded by inserting flat galvanised panels approximately 0.2 m high into the ground to a depth of at least 0.05 m. The plots sizes on the Brigalow Catchment Study were 2.0 m × 1.5 m (3 m²) and those on Belmont Research Station were 2.0 m long by 1.4 m wide (2.8 m²). All plots were on slope gradients less than 5%. A triangular galvanised flume narrowing into a collection pipe was placed downslope to collect runoff water. The difference in plot sizes was due to the width of the flume available at the time of the experiment. Potable water, or rainwater, with low electrical conductivity (<200 mS m⁻¹) was applied using the rainfall simulator described in Loch et al. (2001) at a target intensity of 80 mm h⁻¹ with VeeJet 80100 nozzles. The target intensity was selected to meet a one in 20 year rainfall event as determined by the Bureau of Meteorology Average Recurrence Interval data (Bureau of Meteorology 2022). Runoff was collected in 1-L sample bottles at the commencement of steady state runoff, with water quality samples taken every 3–5 min for 20–30 min. Each runoff sample was timed during collection and weighed after the event to calculate a runoff rate (mm h⁻¹). Water quality samples were analysed at the Department of Environment, Tourism, Science and Innovation Chemistry Centre for total solids and volatile solids (method 2540B) and total suspended solids (method 2540D). Water analysis methods were based on APHA AWWA (2012). Each site was visited once with three rainfall simulation repetitions conducted on each soil type and treatment. A comparison of the observed hydrograph and sediment output per plot and treatment was conducted to validate each plot.

Soil sampling

Surface soil samples (0–0.1 m) were collected pre- and post-rainfall simulation events to determine soil moisture conditions. Soil cores were analysed for air dried moisture content, cation exchange capacity, total organic carbon, bulk density, and particle size fractions of clay, silt, and sand using the methods of Rayment and Lyons (2011). Percent of undispersed particles less than 0.125 mm in the rainfall wet soil surface (P₁₂₅) was determined according to Loch and Rosewell (1992) (Table 1).

Soil erosion modelling: WEPP

The WEPP can be used for multiple-year, long-term simulations as well as single event simulations. For this study, WEPP v2012.8 was run for each experimental plot using the single storm option and plot-specific input files accounting for soil properties, slope geometry (length and slope), ground cover, and rainfall intensity. Soil hydraulic conductivity (K_e) was calibrated in the soil file until the WEPP single storm infiltration rate equalled the observed steady state infiltration rate per plot. The K_e was averaged per soil type and treatment. Steady state was defined as the period when the hydrograph and sediment concentration are constant over time, as distinct from noticeable fluctuations at the beginning of the runoff event. Interrill erodibility in the individual WEPP soil files was adjusted to match modelled sediment loss per unit width (kg m⁻¹) against observations. The observed soil loss was related to observed steady state concentration (g L⁻¹) and converted to kg m⁻¹ via:

where E is the soil loss per unit width (kg m⁻¹), c is the steady state concentration (g L⁻¹), L is slope length (m), and Q is total runoff depth (mm).

Long-term climate data for Belmont Research Station covered a period of 45 years based on Bureau of Meteorology (data from the Rockhampton station. Climate data for the Brigalow Catchment Study were collected onsite, providing a 43-year dataset. These climate data were used to generate the 100-year weather simulation using CLIGEN v5.3.

Soil erosion modelling: RUSLE

The RUSLE equation (Renard et al. 1997) follows:

where A is annual soil loss (t ha⁻¹ year⁻¹); R is the annual rainfall erosivity (MJ mm ha⁻¹ h⁻¹ year⁻¹); K is soil erodibility (t h MJ⁻¹ mm⁻¹); L and S are the slope length and slope steepness factors, collectively referred to as the LS factor; C is the crop cover and management factor; and P is practice factor which includes contouring, terracing, and other conservation practices. The LS, C, and P are unitless, describing the deviation from a USLE unit plot. The USLE unit plot is defined as a length of 22.1 m with a 9% slope, tilled up and down the slope, and retained as bare fallow conditions. Units here are expressed as SI units converted from the original imperial measure by Foster et al. (1981).

For erodibility, the equation for the USLE nomograph K (K_nom) in SI units (t h MJ⁻¹ mm⁻¹) for soils with less than 70% silt used by Rosewell and Edwards (1988) in Australian studies follows:

where M = (percent silt + percent very fine sand) (100 − percent clay), OM is organic matter (%), SS is the soil structure code, and PP is the profile permeability class. The particle size class defined for silt plus very fine sand is 0.002–0.1 mm and for clay is the less than 0.002 mm fraction.

The Loch et al. (1998) nomograph modification (K_m) adjusts the definition of M, where M is replaced by 100(P₁₂₅), P₁₂₅ is particle size distribution less than 0.125 mm, and by using estimated or measured wet sediment density from Loch and Rosewell (1992). Eqn 3 and the K_m modification are used in the GBR Dynamic SedNet and were therefore selected for this study.

In the RUSLE, the ratio of the average soil loss (A) per unit of rainfall erosivity (R) for a unit plot under continuous fallow is represented as K = A/R. In this experiment, LS and C were used to represent a low rill/interrill relationship on a bare plot, and the resulting K-factor follows:

The soil loss A was estimated from WEPP by running 100-year simulations based on the plot-calibrated effective hydraulic conductivity (mm h⁻¹) and interrill erodibility (kg s m⁻⁴) per soil type. The R was obtained from table 2 in Yu (1998) using Station 39083 for Belmont Research Station, and Station 39090 for the Brigalow Catchment Study. A 100-year simulation was run in WEPP to model the annual average soil loss output (A_WEPP) at a slope length of 22 m as reflective of the initial USLE plot at a 9% gradient. The C-factor parameters were adapted from the cover factor equation in Rosewell (1993) where 0.45 is for bare ground (Silburn 2011). The P-factor was set to 1 as it was not applicable here.

Data analysis

RStudio version 4.2.2 was used to analyse and generate graphics of the data using base R, elements of tidyverse including dplyr (Wickham et al. 2023) and ggplot (Wickham 2016) packages. Other packages used to explore and analyse the data include broom (Robinson et al. 2023), multcomp (Hothorn et al. 2023), and car (Fox and Weisberg 2019). Post hoc Bonferroni adjusted multiple comparison analysis was used to identify group differences. Results in Tables 2 and 3 are presented with alphabetical letters where groups not sharing the same letter differ at the 5% level. Performance measures used to compare WEPP models with observations include R², PBIAS, and NSE (Samantaray et al. 2024, 2025). They are defined below, where OBS is observed data, OBS¯ is the mean, and SIM is the WEPP simulated results. The R² is based on the formula used in tidyverse where SS is sum of squares for residuals (res) and totals (tot).

Measured data for the rainfall simulation plots

There were no significant differences for any parameter between natural and artificial scalds on Vertosols and Sodosols (P > 0.7). Based on this result, hydrograph properties for the natural and artificial scalds on these two soil types were combined (Table 2). Minimal to no rilling was observed in the 2-m plots. Runoff started 1.5–4.3 min after the commencement of simulated rainfall, and soils with a higher clay and silt content began runoff earlier. A similar trend was noted for steady state infiltration (mm h⁻¹) where coarser textured soils maintained a higher infiltration rate. Sediment loss peaked early in the runoff events (0–5 min) and gradually declined with time for most plots. Most plots had an observable steady state, with soil loss remaining relatively constant across most times steps except for the Tenosol, which consistently declined for each time step, as did Sodosol₂ to a lesser degree. Steady state observations are summarised for infiltration (mm h⁻¹), sediment concentration (g L⁻¹), sediment generation (g min⁻¹), and runoff rate (mm h⁻¹) in Table 2.

The Dermosol, while higher in clay content (30%) and with a low estimated wet density (1.68 Mg m⁻³), had a high steady state infiltration rate, exceeded only by the sandy loam of the Sodosol₁. The Dermosol and Tenosol runoff depths were the highest at 44 and 43 mm, respectively. Runoff depth results are caveated by the differences in rainfall duration, with most sites receiving 20 min of rain except for the Dermosol and Tenosol with 33 min. All plots receiving ~20 min of rain had runoff depths of 26–32 mm. Runoff initiated earliest at 1.5 min on the Dermosol. The Dermosol had the lowest sediment generation (18.0 g min⁻¹) and concentration (5.4 g L⁻¹) of all soils investigated. The runoff rate from all Dermosol plots averaged 72 mm h⁻¹ with an average rainfall intensity of 84 mm h⁻¹.

The Sodosol₁ was a loamy sand textured soil with a low clay and silt content (9.9% and 3.4%, respectively), and time to runoff from rain initiation was the slowest at 4.0 min. Infiltration was the highest (16 mm h⁻¹) as was sediment concentration (13.1 g L⁻¹). Sodosol₁ had the lowest runoff rate (54 mm h⁻¹). In contrast, Sodosol₂, which had a higher clay and silt content (20.9% and 7.4%, respectively), had the lowest infiltration rate (5.1 mm h⁻¹), the highest runoff rate (83.4 mm h⁻¹), and the lowest runoff depth for plots receiving 20 min of rain (26.3 mm). The Tenosol had the highest sediment rate (51.4 g min⁻¹) and second highest sediment concentration (12.8 g L⁻¹).

The Vertosol₁ showed differences to Vertosol₂, with the higher mean rainfall intensity (91 mm h⁻¹) resulting in a higher runoff rate, depth, sediment concentration, and sediment rate than Vertosol₂. Soil particle size fractions were similar between the two Vertosols; however, Vertosol₂ had higher organic carbon (1.53%) than Vertosol₁ (0.85%). Sediment concentration and sediment generation for Vertosol₂ were low.

All steady state factors of infiltration (mm h⁻¹), runoff rate (mm h⁻¹), sediment concentration (g L⁻¹), and sediment generation rate (g min⁻¹) significantly differed between soil types (P < 0.002). Significant differences were observed between most soil types for runoff depth (mm), primarily influenced by differing rainfall durations and therefore the P-value is not shown here.

Infiltration rates decreased as time increased (Fig. 4), whereas runoff rates increased (Fig. 5). For Sodosol₁ and to a limited degree Vertosol₂, sediment transport increased initially before decreasing as time passed (Fig. 6). Steady state sediment generation was achieved for most soils, with the Tenosol and Sodosol₁ the exceptions. For the Tenosol, and to a lesser degree Sodosol₁, sediment generation was still declining at the end of the rainfall event (Fig. 6).

Fig. 4.

Observed infiltration rate across the runoff events per soil type.

Fig. 5.

Observed runoff rate across the runoff events per soil type.

Fig. 6.

Observed sediment discharge across the runoff events per soil type.

WEPP modelling outputs

Each plot repetition per soil and scald treatment were validated using comparisons of the observed plot hydrograph and sediment output, then combined by treatment and soil type. Artificial scalds were the only treatment observed on the Dermosol and Tenosol. Data modelled in WEPP for infiltration rate (K_e, mm h⁻¹), mean interrill erodibility (K_i, ×10⁶ kg s m⁻⁴), annual soil loss (t ha⁻¹), and the cover factor are shown in Table 3 alongside modelled K_WEPP (erodibility) parameters (t h MJ⁻¹ mm⁻¹). Initial K_WEPP was the result of K = A/R prior to applying Eqn 4.

Total observed soil loss (t ha⁻¹) from each simulated rainfall event was compared to the WEPP modelled soil loss on a plot-by-plot basis to assess WEPP outputs (Fig. 7). The model fit (R² = 0.9, PBIAS = 1.9%, NSE = 0.87) indicates good confidence that the plot calibrations are acceptable.

Fig. 7.

Comparison of observed plot-scale soil loss for each simulated rainfall event against WEPP-modelled soil loss per plot (R² = 0.9, PBIAS = 1.9%, NSE = 0.87, n = 24).

The WEPP mean infiltration (K_e) was calibrated per plot then averaged across each repetition per scald treatment and soil type. Where a plot was not valid, it was dropped from analysis. Hydrological anomalies observed for Vertosol₁ reduced a potential six plots to three, and on Sodosol₁, six plots were reduced to two. The Dermosol included two plots of a potential three. Modelled K_e appeared to follow the same observed steady state infiltration rate trends as field observations, where the sandier textured Sodosol₁ showed considerably higher K_e (20.0 mm h⁻¹). The Tenosol had the lowest K_e of 0.7 mm h⁻¹. The K_e significantly differed between soil types (P < 0.001) (Table 3).

Interrill erosion (K_i) parameters significantly differed between soil types (P < 0.001). Calculated K_i based on Elliot et al. (1989) was tested against modelled WEPP K_i and showed good agreement for most soils and both approaches, with derived K_i not significantly different (P = 0.9). The highest K_i value was for the Tenosol of 2.26 × 10⁶ kg s m⁻⁴ and the lowest was 0.96 × 10⁶ kg s m⁻⁴ for Vertosol₂.

The K_m (modified nomograph) across all soils ranged between 0.032 and 0.048 t h MJ⁻¹ mm⁻¹, indicating moderate erodibility for each soil studied. For K_nom, the range was 0.034–0.040, also indicating moderate erodibility for all soils. When using A_WEPP to substitute into RUSLE, the resultant K_WEPP range for most soils was 0.017–0.028, indicating low erodibility class, with Sodosol₂ the only moderately erodible soil (0.038), which was also the highest K-factor in this study for K_WEPP. For Sodosol₂, K_m is within 8% of the modelled K_WEPP. The Dermosol and Vertosol₁ showed the lowest K_WEPP erodibility (0.017, low erodibility), noting the Dermosol had enough rock cover on the bare surface to reduce the C-factor to 0.39.

The K_WEPP was up to 63% lower than K_m for all soils, except Sodosol₂ where it was 8% higher (Fig. 8). Vertosol₁ and Vertosol₂ were 49% and 58% lower, respectively, when comparing K_WEPP to K_m. Both nomograph-calculated K-factors (K_nom and K_m) indicated moderate to high erodibility for Vertosols (range 0.033–0.042). Sodosol₁ was 17% lower for K_WEPP than K_m. The Tenosol was 42% lower for K_WEPP and the Dermosol 63% lower for K_WEPP.

Fig. 8.

Comparison of K_WEPP and K_m erodibility parameters.

WEPP/RUSLE integration

Limited Australian studies have investigated estimating soil erodibility using WEPP as a substitute for RUSLE, with the exception of developmental work undertaken between 2015 and 2017 as part of the Paddock to Reef Integrated Monitoring and Modelling Program reported in Loch (2017). In Loch (2017), data observed for the Brigalow Catchment Study were modelled using shorter slope length than this study (11 m by 4%) and an extended version of WEPP not publicly available. Loch (2017) used a modified version of WEPP v2012.8, which allows up to 10 different soil size classes as opposed to the three classes used in the publicly available WEPP v2012.8 version, introducing a variability this study cannot account for. No specific methodology or WEPP soil, management, and climate files used to develop K-factors in Loch (2017) are available and were subsequently independently developed for this study.

Internationally, there have been several studies comparing erodibility from the RUSLE and the WEPP. Nearing and Nicks (1998) compared erodibility of 45 soils in the USA and noted that direct comparison between the RUSLE and WEPP erodibilities is complicated by the inclusion of three erodibility-related parameters within WEPP (interrill, rill, and critical shear stress); however, there was no statistical bias between the two model outputs.

Laboratory experiments by Deviren Saygin et al. (2018) investigated the combination of WEPP and RUSLE to estimate the USLE erodibility factor. Deviren Saygin et al. (2018) used laboratory-based rainfall simulation plots to separately calculate interrill and rill parameters to input into WEPP to output an annual soil loss over a 100-year simulation to calculate K. Modelled K in their study was significantly lower than the empirical RUSLE equations for K. The approach by Deviren Saygin et al. (2018) differed to this study by applying a monthly tillage in the WEPP model, where this study considers a permanent consolidated surface as found on grazed hillslopes. That study also used soil transferred to laboratory plots to apply simulated rain, which may not represent the consolidation and organic bonds of in situ soil.

Interrill erosion

As ascertained in Yu and Rosewell (2001), WEPP is a complex process-based model comprising a significant number of initial parameters and sub-routines. Every effort was made to set initial parameters as accurately as possible, based on observed site data. All erosion observed was interrill erosion, as expected for plots with short lengths and low slopes. While this study attempted to isolate interrill erosion specifically as the driver for soil loss, the WEPP model still considers rill formation when scaling the plot data to a 22 m slope length, contributing some rilling to the sediment transport output.

Interrill erosion processes are dominated by raindrop impact and subsequent sediment transport of overland flow (Kinnell 2005; Zhang et al. 2020). When raindrop detachment is greater than the transport capacity of overland flow, the system is transport-limited; if flow is greater than the ability of raindrops to detach soil particles, the system is detachment-limited (Zhang et al. 2020). Both transport- and detachment-limited processes were observed in this study, where initial sediment loss peaked early where raindrop-impact detached soil, declining with time as entrainable sediment supply reduced, and the ability for raindrop impact to disperse particles was less effective as depth of flow increased. In this study, low slopes inhibited flow transport and are representative of the low gradient of grazed lands in the Fitzroy Basin.

Assessing sediment output using WEPP has been tested in several studies. In Australia, Yu and Rosewell (2001) validated soil loss prediction on bare fallow plots, providing confidence in the substitution of A in the RUSLE for soil loss prediction on similarly prepared plots to this study. The WEPP has been successfully applied in various USA conditions to predict sediment yield (Ghidey and Alberts 1996; Ahmadi et al. 2011). Accurately representing the hydrology and climate variables assists in better predictability of soil loss (Laflen et al. 1994), and with known parameters of the rainfall event applied to the plot, measured hydrology, and an appropriately-generated CLIGEN file, the calibration to simulate soil loss was applied with the best available data.

Modelled K-factor

In soils with similar totals of clay and silt (approximately 50% combined total for the Dermosol, Tenosol, and Vertosol₁), observed sediment loss, infiltration, and sediment transport varied significantly. These differences may be influenced by variations in the soil properties of cation exchange capacity and organic matter as captured in the nomograph; however, the variation between the nomograph equations (K_m and K_nom) and the modelled K_WEPP are worthy of discussion. The modified nomograph (K_m), limited to soil properties, indicates that all soils are moderately erodible, with the Tenosol, Dermosol, and Vertosol₁ at the high end of the moderately erodible scale. The K_nom follows a similar trend; however, K_WEPP disagrees with K_m and K_nom.

As an example, the Tenosol had the highest K_m (0.048) and K_nom (0.040) but was in the low erodibility class for K_WEPP (0.028), despite having the highest sediment discharge rate and interrill erosion parameter. Literature shows that Tenosols can vary from low to high erodibility. Silburn (2011) provides a measured K-factor for an Orthic Tenosol of 0.050 in central Queensland’s Nogoa sub-catchment of the Fitzroy Basin, and for a Tenosol in Victoria’s LaTrobe catchment, Vigiak et al. (2011) provides an estimate of 0.015; however, the Tenosol in Victoria had very low clay (1%) and very high organic matter (12%). The Tenosol in Silburn (2011) had 12% clay and low organic matter of 0.6%. The K_WEPP of 0.028 modelled here is a likely representation of in situ, undisturbed soil and seems reasonable for a Tenosol, as is K_m of 0.048.

Vertosol₁ and Vertosol₂ sites were located within the conservatively grazed catchment of the Brigalow Catchment Study. This catchment has a long-term observed soil loss of 0.26 t ha⁻¹ year⁻¹ for the years 1984–2010 (Thornton and Elledge 2021). Substituting K_WEPP (0.018) into the RUSLE equation as K, and appropriate parameters for the catchment (LS = 0.37; R = 2683; C = 0.013), the soil loss estimate was 0.23 t ha⁻¹ year⁻¹. Using the average K_m for Vertosol₁ and Vertosol₂ (0.038), the predicted soil loss output was 0.49 t ha⁻¹ year⁻¹; double the observed data. A K-factor developed by Loch (2017) with a similar WEPP and RUSLE integration based on the same observed dataset (K = 0.023) when substituted into RUSLE, also provided a reasonable estimate for this same catchment (0.30 t ha⁻¹ year⁻¹) compared to the long-term observed soil loss. An analysis of the Modified USLE (MUSLE) by Tiwari et al. (2021) estimated this conservatively grazed catchment had a K-factor of 0.059 which was significantly higher than the WEPP-derived K estimates provided here.

Vertosols are prevalent in the Fitzroy Basin (28% of the grazed land use). Vertosols studied by Freebairn et al. (1996) at the Darling Downs, Queensland, indicate that Vertosols are highly erodible, specifically on cropping systems. The soils in that study had in excess of 50% clay unlike this study. Sediment loss data from the Freebairn et al. (1989) study were evaluated against the USLE, and it was reported that deriving erodibility from the nomograph equation was determined to be unreliable for high clay soils (Freebairn et al. 1989). Measured erodibility from the cropping trial by Freebairn et al. (1989) was used as a point of validation for modelling erodibility using WEPP and RUSLE by Loch (2015b), which showed reasonable agreement, indicating the use of WEPP to derive an annual soil loss average to substitute into RUSLE is feasible.

Aggregate stability is an indicator of susceptibility to erosion (Barthès and Roose 2002; Ben-Hur and Lado 2008), particularly the interaction of clay particles (Lado and Ben-Hur 2004). Vertosols are characterised by high clay content (Isbell and National Committee on Soil and Terrain 2021) and Murphy et al. (2013) makes the point that when wet, Vertosols are highly erodible. At the time of the rainfall simulation study on the Vertosols, the soil moisture prior to applying rain was low (~4%). This study showed variation of sediment discharge between the Vertosols modelled (40.9 and 25.7 g min⁻¹). This may be explained by variations in soil organic carbon (SOC) and clay fractions, with Vertosol₁ exhibiting lower SOC (0.85%) and higher clay content (34.8%) than Vertosol₂ (1.53% and 26%, respectively). Despite these differences, and variations in cation exchange capacity and soil texture, the K_i modelled for the Vertosols was similar (1.06 × 10⁶ and 0.96 × 10⁶ kg s m⁻⁴) and, as a result, so was the K_WEPP simulated under the same climactic conditions (K_WEPP 0.018 and 0.017).

Using the Sodosol₁ dataset, Loch (2017) modelled a 100-year WEPP simulation with a slope length of 11 m by 4% slope with different RUSLE C-factors (0.07 and 0.43) and calculated an average K_WEPP of 0.038, which was higher than the modelled output here (K_WEPP 0.026) by 32%. For Sodosol₂, Loch (2017) derived an average of 0.020 at a LS of 11 m by 4% (C-factor, 0.43). Using a C-factor of 0.45 and a 22 m by 9% LS in this study, Sodosol₂K_WEPP was 0.038, which presents a significant difference. On a larger spatial scale (74 m²), Silburn (2011) calculated a K-factor value of 0.042 for a sandy clay loam, a Brown Sodosol, near Springvale, Queensland, over seven years of observation. This is reasonably comparative to Sodosol₂ but 38% higher than Sodosol₁. The soil texture properties of Sodosol₁ were more similar to the Sodosol studied by Silburn (2011) than Sodosol₂. Local variation is therefore relevant and site-specific K-factors would be useful.

One of the few Australian studies where K-factor was measured were on cultivated modified USLE plots (41 m × 2.44 m, 10% slope, 107 m²) was in Wagga Wagga, New South Wales, where five soils were observed over 4–8 years by Rosewell (1992). That study highlighted the influence of rainfall erosivity variations from year to year on a sandy clay loam (Northcote Principal Profile Form: Dr2.32) similar to Sodosol₂, which generated K-factors of 0.003–0.121. This substantial variability indicates that rainfall erosivity variations over time are relevant, highlighting that long periods of in situ measurement of soil loss are generally advisable to determine the K-factor. In this study, long-term climate variables were captured in the WEPP modelling using CLIGEN, which were based on a long-term rainfall dataset. The variability of the Sodosols across the above studies makes it difficult to assess an appropriate K-factor; however, considering the sodic influence within the Sodosol soil group, erodibility values have the potential to be significant.

This same study by Rosewell (1992) investigated a possible medium clay Dermosol at Inverell, New South Wales (Northcote Principal Profile Form: Uf6.21) on a bare fallow USLE plot 41 m in length at 10% slope, reporting an average K-factor of 0.017. Annual variations observed there were 0.007–0.061. At the same site, a heavy clay Dermosol (Northcote Principal Profile Form: Uf6.21) had an average K of 0.045 (annual variation 0.018–0.070). The Dermosol in this study showed more soil texture similarities to the heavy clay soil at Inverell. A black Dermosol was investigated by Silburn and Bosomworth (2023) in a cropping scenario, and they determined K_i of 3.9 × 10⁶ kg s m⁻⁴, which was significantly higher than K_i estimated in this study (0.59 × 10⁶ kg s m⁻⁴). While both Dermosols had similar clay and silt contents, the loose, unconsolidated condition of the cropping Dermosol resulted in high soil loss under laboratory-simulated rainfall conditions (29.3 t ha⁻¹) compared to in situ plots in this study where the consolidated Dermosol soil averaged a soil loss of 2.4 t ha⁻¹ (Fig. 7). Consolidation may explain the difference in K_i, as well as the cropping Dermosol having considerably higher fine sand content (45%) than the grazing Dermosol of this study (22%). Considering that Dermosols can have high clay or loamy textures, the response to detachment and transport of sediment will be soil and site specific. The K-factor modelled in this study of 0.017 is therefore reasonable and representative of a cohesive surface.

Particle size

The study data suggest cohesiveness, consolidation, and aggregation present in the soil surface influence detachment compared to an erodibility nomograph equation which expects some surface disturbance (e.g. cultivation). Soil aggregate strength, specifically consolidation of in situ soils, influences particle detachment and erodibility, and is integral in erosion resistance (Bryan 2000; Ben-Hur and Lado 2008; Thomaz and Fidalski 2020; Xia et al. 2022; Liu et al. 2023). While particle size is an influential factor, as mentioned by several authors, a variety of other factors increase or decrease resistance to detachment, surface sealing, and entrainment (Alberts et al. 1987; Ben-Hur and Agassi 1997; Ben-Hur and Wakindiki 2004). This was observed here where differing soil property combinations modelled different A_WEPP and K_WEPP values. The Dermosol and Tenosol had high silt percentages of 24% and 31%, respectively, compared to the other soils; however, these soils exhibited vastly different K_i and K_WEPP values. The Dermosol had rock content on the surface and within the soil profile, which may have assisted with infiltration pathways despite a high runoff depth, and thereby inhibited sediment export. Sodosol₁ had low clay and silt (~15% combined) and a higher presence of sand fractions, resulting in a high K_i, whereas Sodosol₂ had a lower K_i yet had the highest K_WEPP of all soils studied. Sodosol₂ had the lowest organic carbon percentage (0.73%), generally an indicator of weaker aggregate strength; however, this alone does not explain the high K_WEPP, as Vertosol₂ and Sodosol₁ had organic carbon less than 1% and modelled low K_WEPP erodibility.

K_WEPP implications

With the exception of Sodosol₂, K_WEPP was consistently lower than K_m and K_nom. This suggests that hillslope RUSLE modelling in Dynamic SedNet, which currently implements K_m, may be overestimating hillslope sediment loss. This has been highlighted in this study where observed and modelled K-factors showed a substantially lower K_WEPP than both erodibility nomograph equations. The nomograph, and indeed the RUSLE, assume a disturbed surface soil not always present on grazed hillslopes when not considering animal traffic itself. Reichert et al. (2001) indicates that increased soil aggregation, consolidation, and surface residue reduce erosion. These conditions are analogous to a grazed hillslope as studied here, so a similar response would be expected.

For the soils considered in this study, which together comprise approximately 70% of the area grazed in the Fitzroy catchment, hillslope erosion could be overestimated, specifically on bare, consolidated ground. With increased ground cover on consolidated soils, the expected soil loss response from hillslope grazing would decrease.

In this study, a storm event was set at one intensity to develop soil loss outputs, then an average long-term erosivity (R) was applied to determine K_WEPP. Undertaking further research over a range of storm event intensities to parameterise WEPP, as done by Silburn and Bosomworth (2023), to develop a range of K-factors against long-term erosivity factors would be of value. The seasonal variability of K, mentioned by several authors (Laflen 1982; Coote et al. 1988; Rosewell 1992; Auerswald et al. 2014; Benavidez et al. 2018), promotes the concept of running repeat simulations with varying rainfall intensities to produce soil loss observations and K-factors more representative of a long-term average. The K_WEPP developed here isolated the dominant interrill erosion process expected on hillslopes of low gradient based on observations from short slope-lengths. With longer slope-lengths and increased rainfall intensities, overland flow and flow power could increase despite the low gradient and introduce the potential for rilling to contribute to erosion. Investigating this process as a subsurface contribution to end-of-catchment would be useful.

While improvements of the K-factor will refine predicted soil loss from hillslopes when applied in catchment sediment modelling, these modifications themselves do not influence observed sediment loads at end-of-catchment. What they do facilitate is an improved understanding of the contribution of the hillslope sediment budget to the GBR lagoon and improved estimation of the spatial distribution of hillslope soil erosion, allowing for focus on either hillslope or other contributing sediment sources at the catchment scale. Further studies, similar to those presented here, will assist in obtaining a wider range of K-factors to capture variances from different soils and rainfall conditions, including the effect of antecedent rain. Where more field-measured K-factors are available for Queensland grazing lands, a new nomograph equation could be derived to better estimate erodibility for GBR catchment soils.