Assessment of Soil Erosion Potential From the Disturbed Surface of Skid Trails in Small Shovel Harvesting System

Forest roads, haul roads, and especially skid trails have been associated with sedimentation and soil erosion risk. Despite the widespread small shovel harvesting system on steep terrains in South Korea, the subsequent risks of deep (rut depth >5 cm) and compact disturbances, and erosion rates in skid trails are largely unknown. Therefore, the main objectives of this study were to compare the soil erosion rate in each skid trail and predict the total soil erosion rate in a small shovel harvesting area. The soil erosion rate was measured at the plot scale (5 × 3 m) in different skid trail parts (bladed skid trail by small-shovel loader passage, BT; and compacted skid trail CT by carrier passage with construction by a small-shovel loader) using a silt fence experiment. In addition, we investigated the applicability of the Water Erosion Prediction Project (WEPP) model to each disturbance. Among all disturbances, the highest erosion rate (average value of 9.13 ± 0.96 kg m−2 4 months−1) was because of CT. The model predictions were over- and under-estimated and showed particularly poor performance where uncovered soil was exposed (less than 1%) to high machine traffic frequency and excavation. Further, the annual soil erosion rates ranged from 11.59 to 28.94 ton ha−1 year−1. The results suggested that the WEPP model could partially validate the soil erosion results, and further research is still required to improve the accuracy of the model.


INTRODUCTION
Forest roads are the most expensive and primary structures for operational access of heavy machines and vehicles necessary for forest management, but these machines can have significant adverse impacts on the environment (Caliskan, 2013;Boston, 2016). They can cause compaction and disturbance of surface and subsurface soil, thereby altering hydrological processes and consequently, causing severe environmental stresses in mountainous terrains because of soil erosion (Fu et al., 2010;Boston, 2016;Kastridis, 2020). Particularly, skid trails, which are considered as temporary pathways for vehicles and which do not have designed drainage or compacted roadbeds, could potentially pose a high erosion risk compared with forest roads after extreme rainfall events (Pierzchala et al., 2014;Safari et al., 2016).
According to the Mountainous District Management Act of South Korea, skid trails have been widely established as a part of forest road networks since 2000 in forestlands because of their low costs and low environmental impacts (Ministry of Government Legislation, 2021). Currently, small shovel loaders and carriers have been commonly used to extract most 2-4 m logs (sawlogs and pulpwood) from the stump to the forest roadside or landing areas (Lee et al., 2019;Lee et al., 2020). During these activities, the skid trails are referred to as vehicle pathways that enable logging machines for transporting logged timber from the harvest sites. Poorly planned forest roads can often cause landslides along road corridors in summer. Moreover, road construction is limited in steep mountain terrains in South Korea. Alternatively, skid trails provide access to vehicles in forest operation sites. They have unpaved road surfaces with widths less than 3 m, in contrast to permanent forest roads, which are designed to have a paved or well-compacted surface. Ji et al. (2020) surveyed skid trails over the entire forest areas of South Korea and indicated that their densities in clearcut harvest units ranged from 76.1 to 184.8 m ha −1 with an average of 108.5 m ha −1 with variations in commercial thinning units from 32.4 to 203.5 m ha −1 with an average of 75.6 m ha −1 .
Heavy equipment traffic can affect subsurface hydrology by increasing the soil bulk density and decreasing the hydraulic conductivity; additionally, it contributes to sediment runoff from unpaved road surfaces (Pierzchala et al., 2014;Safari et al., 2016;Picchio et al., 2020). Machine-induced compaction reduces the soil macropore space and pore connectivity and increases the soil bulk density of the soil surface. Further, eroded soil particles from skid trails eventually reach streams and consequently, degrade the water quality (Clinton and Vose, 2003;Wenger et al., 2018). Particularly, in the surface soils of skid trails at depths of 0-20 cm, small shovel logging systems decreased the soil bulk density by 27-53%, while saturated hydraulic conductivity decreased by 78-98%, compared with natural forest soils . These considerable changes in the soil properties could increase the surface runoff and associated sediment yields from road surfaces during storms (Luce, 2002;Zemke, 2016). Vehicular traffic deforms the trafficked area or causes rutting in pathways (Najafi et al., 2009;Cambi et al., 2015). Ruts may form channel flow along the slope and cause erosion owing to reduced infiltration (Startsev and McNabb 2000;Picchio et al., 2020). In addition, soil erosion can cause the loss of forest canopy and forest-derived leaf litter, reduction of soil organic matter, and decreases interception and infiltration rates (Avwunudiogba and Hudson, 2014;Jourgholami et al., 2019;Picchio et al., 2020). Furthermore, although skid trails occupy only a small part of the total area, they are responsible for a disproportionate amount of total soil erosion from forest areas (Swanson and Dyrness, 1975;McCashion and Rice, 1983;Lewis, 1998;Vinson et al., 2017).
Direct measurements of soil erosion rates are difficult, costly, and time-consuming, but are accurate (Wade et al., 2012;Safari et al., 2016). Since the past few decades, field observations have been conducted to monitor soil erosion in forest road surfaces (Wade et al., 2012;Safari et al., 2016). In contrast, computerbased simulation approaches, ranging from empirical to physical soil erosion models, have been commonly used to estimate soil loss and runoff rates over different geographical regions (Mahmoodabadi and Cerdà, 2013;Gould et al., 2016). Empirical models use location-specific technology and require calibration in diverse areas (Laflen et al., 1991;Raclot and Albergel, 2006). Physical models mathematically describe hydrological and erosional processes and are therefore, applicable to a wide range of soil types, climatic conditions, and land use conditions (Chandramohan et al., 2015). An important and commonly used physical model is the Water Erosion Prediction Project (WEPP) model (Laflen et al., 1991;Amore et al., 2004;Mahmoodabadi and Cerdà, 2013). It has been applied in multiple ecosystems, such as rangelands and forests, at daily, monthly, or annual conditions (Fu et al., 2010;Mahmoodabadi and Cerdà, 2013;Vinson et al., 2017).
It has been successfully applied in various forest conditions, including areas burned by wildfires or intentional fires (Cover et al., 2005;Dun et al., 2009;Elliot, 2013), harvested sites (Elliot et al., 1995;Cover et al., 2005;Haas et al., 2020), forest roads, and skid trails (Elliot and Tysdal, 1999;Fu et al., 2010;Wade et al., 2012). In South Korea, the WEPP was applied to wildfire-burned areas on hillslopes (Park and Shin, 2011) and cultivated uplands (Kang et al., 2004). However, soil erosion is rarely discussed in harvest areas of forests, particularly in small shovel logging operation areas.
Therefore, the main objective of this study was to predict the erosion rate from clearcut harvesting areas using small shovel logging technology. Specifically, the objectives of this study were to: 1) quantify the effects of soil disturbances (undisturbed and shallow disturbance only, US; bladed skid trail by small-shovel loader passes or rutted disturbances, BT; and compacted skid trail i.e., compacted road surfaces, by carrier passages with construction by a small-shovel loader, CT) on sediment yields and 2) evaluate the erosion generation rate for the harvested areas.

Description of the Study Area
The study was conducted in the forests of Joma-myeon, Gimcheon-si, Gyeongsangbuk-do (128°06′48″ E, 36°00′00″ N; Figure 1), South Korea. The climate of this region is characterized by the East Asian monsoon, with a hot rainy summer and a cold dry winter. Since the last decades, the air temperature in the Gimcheon weather station has been varying from 0.1°C (between −5.5 and 6.4°C) in winter to 24.5°C (between 19.8 and 30.4°C) in summer. Two-thirds of annual precipitation occurs in summer from June to August, ranging between 185.6 and 676.2°mm, with an average of 342.5 mm in July, and between 67.1 and 675.7°mm, with an average of 228.4 mm in August.
A 15.1 ha harvest site showed a steep slope gradient ranging from 36 to 47% ( Table 1). The forest type is mixed hardwoodpine forest, and the dominant species are Korean red pine (Pinus densiflora), pitch pine (Pinus rigida), Mongolian oak (Quercus mongolica), and cork oak (Quercus variabilis). The average values of stem diameter at breast height and tree height in the selected study area were 23 cm and 17 m, respectively, and the average stand volume was 142 m 3 ha −1 . The soil was characterized mainly by loamy sand and some sandy loam having 4-11% organic matter content.
The study area was designated for harvest using clearcut harvesting silvicultural treatment from February to April 2017. A motor-manual felling process was used to cut the trees and buck the delimbed trees. Log extraction was performed along the skid trails using a small shovel and carrier, and only merchantable logs were transported to the landing and trafficking areas ( Figure 2). The mass of the small-shovel loader weighed 6.0 ton (including log grapple weight) with a ground contact pressure of 30-40 kPa. The carrier weighed 6.0 ton and had a maximum payload of 6.0 ton of logs with a ground contact pressure of 15-30 kPa. Skid trails were constructed through earthwork by excavating at depths of 30-50 cm and widths of 5 m, and subsequently, compacting the road surface. Further, road density of skid trails was 179.2 m ha −1 . After extraction, the area was visually assessed to characterize the presence and intensity of soil disturbance. Disturbance classes for soil surface disturbance monitoring were described by McMahon (1995) and Curran et al. (2007).

Water Erosion Prediction Project Model Parameterization
The WEPP model was employed to quantify soil erosion generation from the skid trail surface. It is a continuous erosion prediction model that can compute soil erosion along a slope, and sediment yield at a hilltop end, including inter-rill and rill erosion processes (Flanagan et al., 2007). Additionally, the model can represent spatial and temporal distributions of soil  erosion generation in a watershed or on a hillslope. In this study, soil erosion rates were predicted from three disturbed surfaces of skid trails in the timber harvesting area. The major inputs for the WEPP model were climate data, topography data, soil attributes, and management conditions of the study area. Weather data were obtained from a nearby automatic weather station (less than 8 km away from the study area) ( Figure 3). Data comprised daily values for precipitation, time until the peak rainfall intensity, ratio of peak intensity to mean intensity, maximum and minimum temperature, wind velocity and direction, solar radiation, and dew point temperature (Nicks et al., 1995). Other meteorological inputs were generated from the WEPP embedded weather generator (CLIGEN) (Nicks et al., 1995).
A field survey was conducted to obtain the WEPP topographic data for the US, BT, and CT plots. For modeling purposes, hillslope profiles were designated based on the slope length, width, and gradient. The slope length and width were set to 5 and 3 m, respectively, for all plots. The terrain gradients in the US, BT, and CT plots were 30, 48, and 30% of the linear profile, respectively. Further, soil properties, such as sand and clay fractions, organic matter content, bulk density, and initial inter-rill cover, representing the WEPP input data were measured in each plot. Soil bulk density was measured using the core sampling method, and hydraulic conductivity was estimated using a Guelph Permeameter 2800K1 (Eijkelkamp Soil & Water, Giesbeek, Netherlands).
The US plot was considered as the control. Soil and management files were used for each treatment based on the soil type and management practices incorporated in the WEPP, but modified with respect to the above-mentioned soil properties in each trap. The "old forest (20-years forest)" soil and management file was selected for the model run. At the US plot, the sand, silt, and clay fractions varied at 61.7, 7.7, and 4.1%, respectively, and the corresponding values of effective hydraulic conductivity ranged from 37 to 695 mm h −1 .
Further, there were no default land use and management files for representing a bladed trail for small-shovel loader passage in the WEPP database, and thus, "skid trail-every year-disturbed" soil and management file was used with small modifications for the BT conditions. Field observations for sand and clay fractions, organic matter content, and initial inter-rill cover were used as soil property inputs instead of published default values. Sand, clay, and organic matter content were 78.4%-70.0%, 6.5%-60.6%, and 3.0-30.0%, respectively, and factored into the surface cover at 60%. Moreover, the effective hydraulic conductivity varied from 25 to 197 mm h −1 on extremely steep slopes.
"Road-outsloped, unrutted" in the WEPP database was also selected as a guidance for the CT condition. Sand, silt, and clay particle fractions for exposed surfaces in the CT sites were

Silt Fence Measurement and Model Performance Evaluation
After harvesting, sediment traps were immediately installed at three locations (US, BT, and CT) to collect sediment yields. Sediment traps were made of synthetic geotextile fabric with small openings (ranging from 0.3 to 0.8 mm) to hold soil particles (Bugg et al., 2017;Whitman et al., 2019). The traps can retain an average of 83% of sediments and soil eroded from the upslope contributing area (Bugg et al., 2017;Whitman et al., 2019). This technique can be inexpensive data recording method and easy to install (Robichaud and Brown, 2002). However, a major limiting point was the lack of realistic sediment transports and hydraulic flows placed on geotextile fabrics because the amount of suspended sediment is often pass through and unknown (Bugg et al., 2017;Whitman et al., 2019). Therefore, the accumulated erosion rates for a specified period were compared to evaluate the model accuracy with a small-scale study. A 5 × 3 m sediment trap was installed at the centerline of the skid trails ( Figure 4). The sediment traps were operated from May 1 to 31 August 2017, and were periodically weighed to measure the total sediment weight. Subsequently, they were placed into zip lock plastic bags and immediately transported to the laboratory for oven dried weight measuring of suspended sediment.
LR and RE were used to evaluate the accuracy of the WEPP predictions (Moriasi et al., 2007). LR is a straight-line relationship between the observed and simulated data. Moreover, the RE evaluates and compares the accuracy of the simulated data and measured data. Thus, a statistical index was used to determine the agreement between the accuracy of the measured erosion rate and the predicted values.

Statistical Analyses
R statistical package (version 4.0.2) was used for all statistical analyses. Homogeneity was tested using Bartlett's test before comparing erosion rates among three different disturbances (US vs. BT vs. CT). ANOVA and the Games-Howell post-hoc test were used to compare the soil erosion rates among the three different locations. All statistical analyses were performed at a significance level of 0.05.

Effects of Soil Disturbance on Soil Erosion
In the study area, the BT plot (with depths of 5-30 cm) were observed in 29% of the measurement points, and were caused by heavy vehicular traffic. Moreover, the CT plots were observed in 15% of the total observations. The observed sediment yields indicated that both the BT and CT plots produced substantially high erosion rates, whereas the undisturbed surfaces did not influence soil erosion ( Table 2). Soil erosion was not observed in the shallow disturbed areas. BT plots eroded at a rate of 0.55 ± 0.15 (with a range of 0.26-0.71) kg m −2 over a period of 4 months. Further, 60% of the litter cover was disturbed in the bladed skid trails during log extraction and timber transportation. Conversely, CT plots were the primary causes of high erosion rates with an average value of 9.13 ± 0.96 (ranging from 7.62 to 10.90) kg m −2 during the study period. The erosion rates of these plots were up to 41.9 times that of the BT plots because of deforestation and soil disturbances during high traffic intensities for log excavation.
Mean bulk density statistically changed from 1.18 g cm −3 in the US plot to 1.41 and 1.67 g cm −3 in the BT and CT plots, respectively, (p < 0.01) at depths of 0-20 cm ( Figure 5). The average saturated hydraulic conductivity (K sat ) was significantly lower in the BT (56.4 mm h −1 ) and CT (4.7 mm h −1 ) plots than in the US plots (337.5 mm h −1 ; p < 0.01; Figure 6). Organic matter content in the US (6.6%), BT (5.5%), and CT (0.8%) plots showed significant differences (p < 0.01; Figure 7); however, no significant difference was observed between the US and BT plots (p = 0.890). Further, the analysis of variance (ANOVA) test showed that the soil erosion rate was statistically higher in compacted surfaces than in the US and BT plots (p < 0.01). Further, the magnitude of soil erosion in the CT plots was significantly higher than that in the US and BT plots.

Model Accuracy and Soil Erosion Simulation
The WEPP model simulated soil erosion in three different plots (US, BT, and CT) and the corresponding results were compared with the sediment volumes trapped in the silt fences ( Figure 8 and Table 2). Linear relationship (LR) was used to identify the prediction accuracy of the model without calibration. The prediction results indicated that they overestimated soil erosion rates but indicated good model accuracy. However, the relative error (RE) values in the CT plots were satisfactory (11.6-32.1%) based on the information by Moriasi et al. (2007). The predicted erosion rates were considerably higher than the measured values during an increasing trend of the soil erosion rate. Therefore, on comparing the predicted and measured rates, the trends of the WEPP predicted values seemed to be considerably overestimated. This could have possibly occurred owing to the inadequacies in predicting a compact road surface plot only. The RE analysis results indicated that the WEPP model was least suitable in the BT plots because its RE values exceeded the acceptable range. Therefore, the BT plot conditions required calibration to achieve good performance for any disturbed area. Moreover, we revised the default value for the rill erodibility (K r ) parameter of "road-outsloped unrutted" plots. The results indicated that the RE ranged from 6.4 and to 12.8% ( Table 3). The model was more accurate for the trapped sediment data compared with the previously calibrated data. Thus, the WEPP model can be used to predict soil erosion on disturbed surfaces.
The CT surface was a significant contributor to soil erosion during rainfall events. Its density and slope may significantly influence the magnitude of soil erosion (Kastridis, 2020). Therefore, we attempted to analyze the impact of trail density and slope change on the erosion rate under only sandy loam (70% sand and 20% clay) conditions. The CT slope ranged from 30 to 60% and had a CT density of 70-150 m ha −1 . After operation of the small shovel logging system, the BT condition accounted for nearly 28% of the total observations of soil erosion (Lee, 2018). This value was applied to evaluate the soil erosion in the harvesting area.
Our simulation results indicated that annual soil erosion rates potentially varied from 11.59 ton ha −1 year −1 for low CT density (70 m ha −1 ) and a slope of 30% to 28.94 ton ha-1 year-1 for high CT density (150 m ha −1 ) and a slope of 60% during an annual precipitation of 1,020 mm and maximum rainfall intensity of up to 33 mm h −1 after the small shovel logging system (Figure 9). This implied that the increase in the trail slope and density increased the soil erosion generation. Thus, the trail slope and density significantly contributed to the total erosion rates, particularly in typical small shovel harvesting systems.

DICUSSION
During forest management activities, most sedimentation and soil erosion events may occur severely because of inappropriately developed forest roads and skid trails. Further, machine traffic can alter the physical and hydrological properties of hillslope soils by increasing the bulk density and decreasing infiltration. Studies on soil erosion rates associated with clearcut operations in Korean forests have not been conducted previously. This study attempted to determine the effect of skid trails associated with small shovel harvesting systems on the soil erosion generation rate and validate the WEPP model at the plot scale. When soils were not covered by vegetation and exposed to high machine traffic frequency, the average soil erosion rates were approximately 20 times higher than the bladed skid trail formed by the passage of small shovel loaders. The WEPP model showed the highest accuracy for the BT plots based on the LR and RE, while it did not give satisfactorily accurate results for the CT plots. In the CT plots, the modified model with Kr outperformed the default predictions. After model calibration, the annual soil erosion rate was estimated to range from 11.59 to 28.94 ton ha −1 year −1 under various conditions. Our results provided significant information for  improving the implementation of forest management practices associated with small shovel harvesting operations.
The disturbances of skid trails may be restricted to reducing soil erosion in comparison to undisturbed areas. Elliot (2013) noted that an artificial linear structure (i.e., forest roads and skid trails) can cause a perpetual erosion problem (0.10-10.00 kg m −2 ), whereas undisturbed forests cause low erosion (<0.01 kg m −2 ). The soil erosion of skid trails constructed using a bulldozer (John Deere 450E; 12.39 kg m −2 ) was considerably higher than mulch treatments (0.27 kg m −2 ; Wade et al., 2012). Similar results were observed in this study, with the highest erosion rate in the CT plot. In contrast, Vinson et al. (2017) reported that skid trails, which are defined as vehicle traffic without any trail construction, eroded at a rate of 1.36 kg m −2 Safari et al. (2016) observed that the erosion rates in wheel tracks and central areas between tracks were 0.30 and 0.26 kg m −2 , respectively, with a significant difference in the erosion rate between the plots (p < 0.001). Further, soil erosion in the BT plot showed a similar trend to that reported in another skid trail study (Safari et al., 2016).
The high level of compaction due to high traffic intensities can be explained by the high soil erosion rates (Zemke, 2016). Additionally, Zemke (2016) reported that K sat is an important factor that directly affects infiltration and soil erosion rates. Moreover, erosion has been observed to be caused by surface conditions, such as fortified and unfortified surfaces (Jordán and Martínez-Zavala, 2008;Zemke, 2016). The resulting removal of ground cover and compaction increases the risk to erosion (Notcliff et al., 1990;Parsakhoo et al., 2009;Elliot, 2013;Wenger et al., 2018). Furthermore, according to Kastridis (2020), unpaved roads increase the erosion risk potential ability than undisturbed forests, grasslands, and agricultural lands. Thus, the CT plot may extensively influence soil erosion in timber harvesting areas.
WEPP is a physical model and shows more accuracy than an empirical model; additionally, it may overestimate soil erosion rates. These findings were consistent with many previous studies, such as those by Covert et al. (2005), Dun et al. (2009), Wade et al. (2012, Vinson et al. (2017), and Wu et al. (2018), who examined the WEPP performance. The findings of these models indicated that model simulations could not only potentially overestimate deep percolation but also underestimate subsurface lateral flow in forests due to complex root systems, shallow depth soil layers, and steep rocky slopes. Our results showed that the WEPP prediction error rate increased linearly with increasing magnitudes of erosion in compacted road surfaces. Similar results were obtained by Mahmoodabadi and Cerdà (2013). The inadequate performance can be primarily due to road surfacing material, such as road building and soil rock content, and absence of management practice files for compacted road surfaces (Vinson et al., 2017). In addition, this model overestimated the rill erosion for a short (up to 10 m) hillslope, with high rill densities, than long hillslopes (Wu et al., 2018). Therefore, for further applications of the WEPP model, model predictions need to be calibrated and validated with locally obtained data.
Model calibration is necessary to improve the accuracy of WEPP model predictions. According to Raclot and Albergel (2006), Foltz et al. (2011), and Mahmoodabadi and Cerdà (2013), the WEPP model has been validated well, whereas the erosion parameter sets require specific location data for improving the accuracy. For example, Dun et al. (2009) and Mahmoodabadi and Cerdà (2013) emphasized that the model requires further modifications to properly perform hydrological processes for forest applications. During simulation modeling for a hillslope area, rill erodibility could be a primary influential indicator than inter-rill erodibility and critical shear stress (Nearing et al., 1990;Mirzaee et al., 2017). Our results indicated that after Kr value calibration, the WEPP model showed better performance for the observed values than default prediction values, only in the CT plot. The Kr value for the skid trails showed a wide range associated with vehicle traffic after construction in harvested area in the Idaho State, US (0.0003-0.0005 s m −1 (Covert et al., 2005) and 0.0076-0.0005 s m −1 (Foltz et al., 2008)). In this study, K r was modified from its default value (0.0004 s m −1 ) to 0.0003 s m −1 because the model tended to overestimate. During the application Frontiers in Environmental Science | www.frontiersin.org June 2022 | Volume 10 | Article 756848 8 of the WEPP model, the predicted erosion rate, associated with a skid trail (compacted road surfaces) formed by carrier passages with construction by a small-shovel loader, may be perilously overestimated. Further, the primary reason for the low model accuracy can be the lack of management files for the CT plots. Therefore, WEPP predictions may require specific location data associated with soil and management files to improve accuracy.

CONCLUSION
In conclusions, this study was the first attempt of its kind to understand the soil erosion generation rate associated with small shovel harvesting systems in different parts of skid trails, such as undisturbed and shallow disturbance only, bladed skid trails by small-shovel loader passes, and compacted skid trails by carrier passages with construction by a small-shovel loader. After designing nine experimental plots for 4 months using the WEPP model, we analyzed the model accuracy. The magnitude of soil erosion rate in skid trails was significantly higher within uncovered surface conditions caused by excavation and high traffic intensities than in the US and BT plots. Based on the LR and RE analyses, the WEPP model proved to be insufficient for use in CT conditions, because it overestimated all values. After calibration with rill erodibility in only CT, values were closer to the observed value than the default predictions. Although the validation of the WEPP model at different locations has not been tested extensively, our results may be well supported by previous studies for further applications of the model. Future research is required to monitor the sedimentation and soil erosion rate in the study area and evaluate the WEPP model calibration and validation under different conditions, such as slope length, soil rock content, soil texture (i.e., loam, sandy loam, and clay), and moisture content.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
EL, SE, and QL contributed to conception and design of the study. EL and SE organized the database. EL and ES performed the statistical analysis. EL wrote the first draft of the manuscript. E, SE, and QL wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
This research received no external funding.