Impacts of wood extraction on soil: assessing rutting and soil compaction caused by skidding and forwarding by means of traditional and innovative methods

Intensive forestry operations may cause soil compaction, plastic soil disturbances and rutting, which are responsible for undesirable effects on soils, vegetation and water bodies. Despite the numerous studies aimed to identify the main factors affecting soil damages, it still remains unclear whether wood extraction methods and driving direction (uphill or downhill) may affect the impacts of forest machines. This research analyses soil compaction and soil penetration resistance as well as rutting from forwarding and skidding using the same farm tractor in up- and downhill wood extraction. Rutting was estimated by 3D soil reconstruction derived by portable laser scanning (PLS) and close-range photogrammetry using structure for motion (SfM). Our findings showed that the direction of extraction did not affect soil damage severity during forwarding on a 25% slope. On the contrary, in order to reduce soil compaction, downhill skidding is preferable to uphill skidding. The results showed that the pressure on the ground caused by vehicles can be distributed horizontally, thus affecting also the soil between the wheel tracks. The soil bulk density inside the tracks after 10 forwarding passes increased by 40% and with 23% between the wheel tracks. The soil displacement in skidding trails (7.36 m3 per 100 m of trail) was significantly higher than in forwarding (1.68 m3 per 100 m of trail). The rutting estimation showed no significant difference between the PLS and SfM methods, even comparing the two digital surface models (DSMs) obtained, even if photogrammetry was preferred for technical and practical reasons.


Introduction
Soil disturbance is an unavoidable consequence of timber logging, but the severity of its impact is variable and can be managed through good planning and practices (Ares et al. 2005). According to the SFO (sustainable forest operation) concept, the minimization of soil impacts is a key task to improve the environmental efficiency of timber logging (Marchi et al. 2018). Intensive forestry operations may cause soil compaction, plastic soil disturbances and rutting, which are responsible for undesirable effects on soils, vegetation and water bodies (Cambi et al. 2015). Ground-based logging operations can negatively affect soil physical characteristics, reducing porosity while also increasing bulk density and resistance to mechanical penetration (Siegel-Issem et al. 2005;D'Acqui et al. 2020;Lee et al. 2020). Moreover, soil compaction indirectly influences tree growth and regeneration due to both physical root damage and reduced soil permeability (Cambi et al. 2018a(Cambi et al. , 2018bJansson and Wasterlund 1999;Mariotti et al. 2020;Sirén et al. 2013;Solgi et al. 2019;Sugai et al. 2020), which may lead to deficiencies of oxygen, water and/or nutrients (Batey 2009;Lee et al. 2020) with recovery processes that may take several decades (Bottinelli et al. 2014;Jourgholami et al. 2020). Rutting and other soil disturbances can also disperse pathogenic fungi, alter microbiological processes (Frey et al. 2009;De Wit et al. 2014;Cambi et al. 2017) and mobilize heavy metals due to the increase in surface water flow (Frey et al. 2009;De Wit et al. 2014;Eklöf et al. 2014). Depending on logging conditions (e.g. soil condition and type, wood extraction method, machine characteristics and operator skills), the surface affected by disturbance within the logging area may range widely from 10 to 87% (Spinelli et al. 2010;Marchi et al. 2014;Naghdi et al. 2015).
Several studies were published on the impact of logging on soil in recent years; in this field, the most important factors affecting compaction and rutting on forest soils and their relationships and interactions, such as slope gradient, driving direction (uphill vs. downhill), extraction method (skidding vs. forwarding), still need to be examined in depth. This is fundamental to be able to plan effective mitigation strategies. Some studies highlighted the effect of slope on soil compaction, finding that the higher the terrain slope, the higher the impact on soil (Jamshidi et al. 2008;Jourgholami et al. 2014;Naghdi et al. 2020). Previous studies have also investigated the impact caused by skidding and forwarding, highlighting that the extent of surface disturbance was higher after skidding, whereas compaction showed contrasting results between the two systems (Lanford and Stokes 1995;Deconchat 2001;Gondard et al. 2003;Cambi et al. 2018a). Compared to forwarders, skidders are driven in a denser extraction trail network (Han et al. 2009) with lower loads, thus explaining the greater extent of surface damages. In forwarding, the impacts on soil are related only to the contact of tyres on the ground (rut formation), whereas in skidding, the semi-dragged logs also contribute to soil impact (soil displacement). In general, impacts on soil due to skidding have been more closely studied than forwarding. In addition, neither data nor studies have been found that compare the effects on soil caused by skidding and forwarding considering a similar operational context, which highlights a gap in the literature.
In recent years, innovative tools that are useful in assessing soil surface disturbances have become available. For the measurement of rutting which relied on classical manual measurements (Jester and Klik 2005) are now available some innovative techniques such as three-dimensional ground reconstruction derived by photogrammetry (Pierzchała et al. 2014(Pierzchała et al. , 2016Haas et al. 2016;Marra et al. 2018;Talbot et al. 2018) or laser scanning (Koreň et al. 2015). In particular, portable laser scanners (PLS) have started to be used in forestry for rutting estimation (Giannetti et al. 2017). However, sizeable knowledge gaps are evident concerning the efficacy of these newer methods (i.e. PLS and close-range photogrammetry) in the impact assessment of forest operations.
The objective of this research was to study the effects of timber extraction on steep terrain by measuring bulk density (BD), penetration resistance (PR) and surface disturbances such as rutting, considering both the extraction method (skidding and forwarding) and driving direction (uphill and downhill). Also, the application of two new methods such as 3D soil reconstructions derived by both portable laser scanner and close-range photogrammetry using Structure for Motion (SfM) were tested for measuring rutting and to identify the most suitable for applications in a forest context.

Site description and experimental design
The fieldworks for this study were carried out from June 15 to July 5, 2018, in the Rincine forest, a public forest property located in the north-east part of Florence province, 40 km outside the Florence town (central Italy, N 43°52′; E 11°34′; 400 m above sea level). Considering climate data from the last 30 years, the climate of the study site has been classified as Mediterranean (Köppen-Geiger classification), characterized by a hot, dry summer (with January as the coldest month and August the hottest). In the period 2017-2018, the local mean annual temperature was 9.2 °C (1.5 °C in the coldest months and 17.8 °C in the warmest), and the mean annual precipitation was 924 mm, with the maximum 1 3 in November and the minimum in July. The soil of the study site developed on Lower Miocene-Oligocene sandstone and was classified as Dystric Cambisol based on the World Reference Base for Soil Resources (IUSS Working Group 2014). The study site was identified within a 50-year-old high stand of Douglas Fir, located 900 m above sea level, characterised by a density of 800 stems·ha −1 and trees with 40 cm of average diameter and 25 m of average height. At the time of the study, no logging activities had been conducted within the study area in the last 40 years, thereby avoiding any influence of previous logging operations. Moreover, no silvicultural treatments were made during the study to avoid additional and unpredictable variables, focusing only on the extraction operations. These were made using a set of Douglas Fir logs extracted in a similar parcel next to the study area. Four trails were designed on the ground in the direction of the slope. Two trails were used for forwarding and two for skidding; in both extraction methods, one trail was used uphill and the other downhill ( Fig. 1). All trails had the same slope (25%). On these trails, wood extraction cycles were simulated using a four-wheel driven agricultural tractor adapted to work in forestry for skidding (equipped with a winch) or forwarding (equipped with forest trailer) Douglas Fir logs. These machines are very common in forest operations in Italy, especially in peninsular forests, and more generally throughout the Mediterranean area. The tractor was a New Holland T6050 (power 93 kW), the winch a Farmi JL61 and the trailer a Zaccaria ZAM 140 Forestal TC Super. The trailer had two axles with two wheels each; the front axle was powered by the mechanical transmission from the power take off (PTO) from the tractor, while the rear one was a not powered self-steering axle. Technical characteristics of the machines are summarized in Table 1. The contact areas shown in the table were measured by means of a rope pulled tightly around the portion of the tyre on the ground and assuming a circular contact patch (Cambi et al. 2016). Ground contact pressure was calculated as a ratio between the machine mass and the contact area (Table 1).
The tractor (in both configurations, with trailer or winch) moved on the four new trails designed within the stand to monitor physical parameters and rutting before and after machine traffic. The tractor operator was an experienced forest worker who drove as in a standard operation along the trail. Three rectangular plots (length 3 m, width 6 m) were selected and marked on the ground along each straight trail section (Fig. 1).
Volume and number of logs loaded in skidding and forwarding are reported in Table 2. The volume of each log was calculated applying Cavalieri's formula: where: V l = volume of the log; S t = section at top-end; S b = section at bottom-end; S 0.5 = section at half length; h = length. Regarding the weight of logs, it has been In skidding, the largest ends of the logs (bottom-end) were closest to the tractor and suspended against the base plate of the winch at an average of 0.8 m above the ground, while the smallest end (top-end) was dragged on the ground. Considering that the load capacity in forwarding was approximately 4.5 times higher than the load capacity in skidding, the comparison of skidding and forwarding was made by taking into consideration both the number of passes and the total volume moved. In detail, during the experimental study, the tractor equipped with a loaded trailer (forwarding) passed 10 times on both the trail driven uphill and the trail driven downhill. The tractor equipped with a winch passed 10 times with logs on each trail and, after an intermediate data collection, made 35 more passes (45 passes in total). On each trail, all passes were made with load. In this way, the same volume of wood was moved with the two systems (≃127 m 3 ), and it was therefore possible to compare the impacts related to both the same number of passes and the same volume of wood transported.

Physical soil parameters
Impact on soil physical parameters and rutting were measured on all plots of all trails, before and after 10 and 45 (only for skidding) passes, applying different methodologies. The physical parameters measured for determining soil compaction were bulk density (BD) and penetration resistance (PR).
Before machine traffic, three soil samples were collected approximately one meter from each identified trail, on both the left and right sides, for a total of six samples per plot.
After the machine passes (10 forwarding passes and 10 and 45 for skidding), three soil samples were collected inside both left and right wheel tracks (InR) and three samples were collected between the tracks (BtwR), for a total of nine samples per plot (Fig. 2). All soil samples were collected from the top 10 cm of mineral soil layer, using a metal cylinder (7.5 cm inner diameter and 10 cm height). BD was calculated as the dry-weight (after a treatment at 105 °C per 48 h) of the soil sample divided by the volume of the sampler (Picchio et al. 2009). In addition, six soil samples were collected daily from the trails (outside the plots) in order to monitor the soil moisture over time (soil moisture content on a dry-weight basis). The PR was measured by a hand cone penetrometer (Field Scout SC 900). A cone with a diameter of 12.7 mm and cone angle of 30° was used to register the average penetration resistance in the first 10 cm of the soil (at each 2.5 cm of soil depth). In each plot, PR was measured following the same protocol as BD (3 samples each in InR and 3 in BtwR) as an average of three values collected in each sample point.

Rut measurements
Two methods to implement a 3D soil reconstruction, derived by both portable laser scanner (hereafter PLS) and  close-range photogrammetry using structure for motion (hereafter SfM), were applied to measure rutting caused by machine traffic on forest soil. Measurements and data acquisition took place in the same plots used for collecting the other soil physical parameter data. Each method may provide accurate information about changes to the ground surface caused by the passage of forest machines by generating 3D representations of the trail section before and after the machine passes. The 3D soil reconstructions derived by both PLS and SfM were used for determining the volume of ruts caused by machine wheels and/or by dragged logs after skidding and forwarding, thereby allowing for the calculation of differences between digging and carryover soil volume.

Data collection
The application of SfM is based on a series of images collected in the field. In our study, these images were collected using an RGB reflex camera: a Canon EOS 1300D model with 18 MP of image resolution and 5184 × 3456 pixels. The camera, which has a focal length of 18 mm, was mounted on a tripod 3 m in height with an angle of 45°. Image acquisition followed a specific scheme, previously described in Marra et al. (2018), designed to collect several images every half meter in different directions. All the images collected were checked directly in the field, immediately after collection, by visual interpretation to detect eventual problems related to light, saturation and blurriness, recollecting lowquality. Before the acquisition of the images, nine ground control points (GCPs) were positioned on each plot for image geolocation (Fig. 3). GCPs were composed of Schneider's coded target (Schneider 1991), an automatically recognizable feature generated with Agisoft Photoscan® software (Agisoft 2016). In detail, each GCP was a 144 cm 2 square aluminium plate with a specific printed image. The image was characterized by a central dot (radius of 15 mm) surrounded by a code band (12-bit pattern) with bit positions at equally spaced angular intervals (black or white). A local system of coordinates was created and the X, Y, and Z coordinates of each GCP were measured and used to align the models acquired at different phases.

Pre-elaboration
The SfM technique was applied to photo analyses to obtain 3D georeferenced point clouds. The analysis has been implemented in several steps using the Agisoft Photoscan® photogrammetric software package: (i) image import, (ii) image alignment, (iii) georeferencing, (iv) optimization of image alignment and (v) creation of the point and dense clouds. After a rough alignment of loaded photos, the camera positions were optimised through the automatic detection and matching of GCPs using a specific functionality of Agisoft Photoscan®. This guaranteed sub-pixel accuracy without the need for human interventions. Thanks to camera optimisation and GCP detection, the single dense point Scheme of soil sample and data collection on the trails for physical soil characteristics after machine passes. On each trail, the data were collected within the track (InR) and between the tracks (BtwR). The circles indicate soil samples collected for bulk density; the diamonds indicate penetration resistance measurements Fig. 3 Position of the Ground Control Points during laser scanning (above) and close-range photogrammetry (below) data collection cloud was produced to determine the 3D soil reconstruction of each plot point at each collection time. All workflow was elaborated using Python 3.5 (Van Rossum and Drake 2020) as a scripting engine and supported by Agisoft Photoscan®.

Data collection
A lightweight, portable ZEB1 HMLS consisting of a 2D laser scanner (Geoslam: Ruddington, Nottinghamshire 2016), combined with an inertial measurement unit (IMU), was used to scan the plots within the study area collecting spatial data (acquisition speed 43,200 points/sec). The reported outdoor operative laser range is 15-20 m around the instrument (Bosse et al. 2012;Giannetti et al. 2017) with a scan ranging noise of ± 30 mm. In each plot, four GCPs were mounted as spherical targets (diameter = 0.14 m) at the top of 1.5 m poles (Fig. 3) fixed in the ground at the corners of the plot. The X, Y and Z coordinates of each GCP were measured in the same systems of coordinates implemented for SfM, thereby guaranteeing the correct overlapping of results for comparison. To collect data in the plots (Fig. 4), the operator walked slowly (approximately 30 cm s -1 ) holding and oscillating the instrument at breast height (1.40 m above ground). The route inside the plot area was covered by walking the entire plot surface back and forth along straight lines (spaced 0.5 m apart), beginning and ending the survey at the same point (Ryding et al. 2015). In order to avoid shadow zones and to obtain the highest scan density of GCPs, special care was taken to fully scan the spherical targets.

Pre-elaboration
The 3D georeferenced point cloud was calculated from raw data by an online service (Geoslam: Ruddington, Nottinghamshire 2016) provided by the producer of the ZEB1 laser scanner. This procedure uses a novel 3D simultaneous localization and mapping algorithm (SLAM) to combine the 2D laser scan data with the IMU data to generate accurate 3D point clouds (Bauwens et al. 2016). The resulting point clouds were loaded in CloudCompare (CloudCompare Version 2.11 2020) in.las file format and then 'cleaned' to detect the plot and remove non-soil objects, such as stumps and trees. This process included the following steps: (i) identifying the plot in the point cloud by visual interpretation of GCPs, (ii) manually tracking and segmenting the plot, and (iii) georeferenced processing of each GCP.

Co-registration of point clouds
The point clouds for both methods used were elaborated by CloudCompare software to obtain the digital surface models (DSM) before and after forest machines passes in all plots. First, the soil surface conditions before (i.e. reference cloud) and after (i.e. compared clouds) machine traffic were coregistered in CloudCompare using the GCPs and analysed using the 'point pair-based alignment tool'. A fast control was conducted to assess the reliability of the automatic cloud's overlap. In addition, minimum errors in positioning were eliminated through a calibration of the coordinate system considering the root-mean-square error (approximately 4 cm) for vertical and horizontal displacement. Finally, the point clouds were rasterized considering the average Z coordinate to create high-resolution (pixel = 1 cm) DSMs for each plot.

Rut measurement estimations
The differences between the DSMs before and after machine traffic were implemented in CloudCompare in order to determine the volumes of ruts (volume reduction = VR) and bulges (volume increased = VI) caused by machine wheels or dragged logs. In addition, all the DSMs derived by SfM and PLS had the same resolution (601 × 301 pixels) to compare accuracy, pixel by pixel, with the R-cran raster packages (Hijmans et al. 2020).

Physical parameters and rutting: analysis of data
Statistical analyses were performed using R software (R Development Core Team 2020). After checking for normality (Pearson chi-square test) and homogeneity of variance (Bartlett test), multi-way ANOVA was applied to the BD in order to test the soil compaction effects by machine type, impacted zone (InR or BtwR), driving direction (uphill or downhill), number of passes and total volume moved. A post Fig. 4 Operator during laser scanner survey, holding and oscillating the instrument at breast height (1.30 m above ground) while walking hoc LSD test (least significant difference test) was applied to rank the results between dependent significant variables. Statistical differences in soil moisture at the time of each field day were tested by means of a t test. The Kruskal-Wallis nonparametric multiple-comparison test (Dunn's test) was used to PR data because data distribution was not normal and variances were not homogeneous. Regarding the analysis on rutting, to verify the accuracy of DSMs obtained with SfM and PLS, regression analyses and the root-mean-square error (RMSE) were performed. Finally, the soil volume changes (VI and VR) were compared by applying a oneway ANOVA in order to understand the effects of the wood extraction methods and to compare the volume estimation with SfM and PLS.

Physical soil parameters
In the study site, the average soil moisture was 18% throughout the entire study period. BD and PR were significantly affected by machine passes. During forwarding, a larger load than skidding was transported in each pass. However, a higher number of skidding passes was required to extract the same wood volume than in forwarding (10 passes in forwarding and 45 in skidding, to obtain 127 m 3 wood volume moved). In detail, the weight of the vehicles and their loads for each pass was 17,962 kg in total for forwarding (7022 kg of load) and 8294 kg for skidding (1549 kg of load). The total weight transported on the ground of the study site after the woods extraction, was higher for skidding (373,246 kg) than for forwarding (179,624 kg) on the study site.

Bulk density-BD
The BD measured within the plots (both InR and BtwR) was significantly higher than in the undisturbed soil before the trials (Table 3).
Comparison of BD considering the same number of passes in skidding and forwarding After 10 passes, in both skidding (28 m 3 of wood volume transported on the trails) and forwarding (127 m 3 of wood volume transported), a significant increase of the BD values was found (Table 3). BD was not affected by the driving direction; however, it was significantly affected by the logging method applied and the impacted zone (InR or BtwR). In general, skidding showed significantly lower values of BD in InR than forwarding, while similar BD values were recorded in BtwR. Finally, when forwarding moving both uphill and downhill, the increase in BD InR (40% higher than control on average; range 25%-52%) was significantly higher than the increase in BD BtwR (23% higher than control on average). After skidding operations, an average increase of 30% (range 17%-40%) was recorded in both BD InR and BtwR.
Comparison of BD considering the same wood volume moved by both skidding and forwarding Increasing the number of skidding passes to simulate extraction of the same wood volume by forwarding (127 m 3 ) affected the measured values of BD. Comparing the effects of the two extraction methods, the results showed similar BD values when operating uphill, while when working downhill, skidding showed significantly lower values of BD in InR-and significantly higher values of BD in BtwR-than forwarding (Fig. 5). Thus, the highest average values of BD were obtained in InR, in both uphill and downhill directions for forwarding, but only in uphill for skidding (approximately 40% higher than undisturbed soil). On the other hand, the effects on soil in BtwR were lower than in InR, except when using the tractor with winch downhill, in which case the dragged logs caused a greater increase of BD in BtwR than the tyres in InR. As reported in Table 4, the analysis of variance showed a strong influence of both the wood extraction methods and the impacted zone (InR/BtwR) on soil compaction after machine traffic, in terms of BD. Moreover, a significant effect of driving direction (uphill vs. downhill) was recorded only after 45 passes. The statistical difference between InR and BtwR was strongly significant both after 10 passes and 127 m 3 moved. Analysing multiple interactions for BD, significant ones were found for wood extraction method and driving direction and between driving direction and impacted zone only when considering the same wood volume moved, while those interactions were not-significant when considering the same number of passes. The interaction between wood extraction method and impacted zone (BtwR vs. InR) was stronger when considering the same number of passes than the same volume moved. Finally, the interaction among wood extraction method, driving direction and impacted zone was strongly significant when considering the same wood volume moved.

Penetration resistance-PR
Comparison of PR considering the same number of passes in skidding and forwarding After 10 passes of tractor with trailer (i.e. 127 m 3 moved), only the values of PR recorded within the tracks (InR) were significantly higher than the control (untrafficked) ( Table 5). The PR values between the tracks (BtwR) did not change significantly. When considering the forwarding direction, significant differences were found between InR and BtwR when the tractor with trailer was driven uphill. Any statistical difference was recorded in skidding after 10 passes (i.e. 28 m 3 moved) for both driving direction and within (InR) or between tracks (BtwR). The comparison of skidding and forwarding did not show any significant difference, neither for driving direction nor InR/ BtwR. The highest value of PR was recorded in forwarding uphill InR. This value was significantly higher than that recorded in skidding downhill InR and BtwR.
Comparison of PR considering the same wood volume moved by both skidding and forwarding Further increments in PR (both InR and BtwR) were recorded in skidding trails after 45 passes (i.e. 127 m 3 -the same wood volume moved by forwarding). Increasing the number of passes, PR did not show any significant differences among skidding treatments (Fig. 6). Comparing skidding and forwardingand considering the same wood volume moved-the results highlighted that the lower PR values were recorded in for- warding BtwR, both in the uphill and downhill directions.
The highest values were recorded in InR after forwarding and in all the skidding treatments. The PR value of forwarding BtwR uphill was significantly lower than skidding BtwR, both in the downhill and uphill directions (Fig. 6).

Rut measurements
The application of two methods-SfM and PLS-to measure the ruts due to wood extraction produced two different 3D reconstructions for each plot implemented in DSMs.

SfM dense cloud
The 3D reconstructions derived by SfM close-range photogrammetry had approximately 2,600,000 points for each plot (144,400 points/m 2 ). The errors of co-registration were always less than 0.13 mm for both the x and y coordinates, and less than 0.31 mm in the z coordinate thus negligible for our purpose.

PLS point cloud
The 3D reconstruction derived by PLS had approximately 1,600,000 points for each plot (88,800 points/m 2 ). The errors of co-registration were always less than 0.21 mm for both the x and y coordinates, and less than 0.31 mm in the z coordinate thus negligible for our purpose. A total of 10 plots were elaborated with PLS instead of the 12 that had been planned, because two files containing raw data (one forwarding downhill and one skidding uphill-both after 10 passes) revealed corruption during the post-processing phase.

DSMs: SfM vs. PLS
The correspondence between the SfM and PLS methods can be easily shown in the pixel-by-pixel comparison of the resulting DSMs obtained by overlapping them. The overlap of each plot showed an RMSE ranging between 2 and 4 cm; in the 97% of the plot area the difference between the two methodologies was less than 2 cm. A spatial representation of this comparison is reported in Fig. 7. In addition, this comparison indicates a significant relationship between the 3D reconstruction data of the two methods (R 2 range between 0.97 and 0.90; p value < 0.000). Moreover, the comparison between these two methodologies focus to measure ruts did not show a significant difference (p value < 0.001), which reveals a similar quality and reliability of results between both methods with a good correlation (R 2 0.80). Examples of 3D soil reconstructions derived by PLS and SfM are shown in Fig. 8.

Quantification of rutting from DSMs
Because no significant difference was recorded in the estimation of rutting using photogrammetry and PLS (Fig. 9), only the results deriving from photogrammetry are reported below. However, significant differences in VI were found in relation to wood extraction method when considering the same wood volume moved. In detail, VI after skidding (7.36 m 3 /100 m of trail) was almost four times higher than after forwarding (1.68 m 3 /100 m of trail). Conversely, the VR derived by skidding was −8.88 m 3 /100 m of trail on average, which was similar to VR after forwarding (−8.51 m 3 /100 m of trail). No significant differences were found in VR and VI considering the driving direction in both extraction methods.
A strong and significant positive relationship between VR and soil compaction was found (R 2 0.704; p value < 0.001) after 127 m 3 of wood volume moved: the higher the VR, the higher the soil compaction. In contrast, the relationship between PR and BD showed a low coefficient of determination (R 2 0.20; p level < 0.000).

Discussion
In this study, we compared the impacts (compaction and rutting) on forest soil caused by skidding and forwarding, considering both the same number of passes by loaded machine (10 passes for both skidding and forwarding) or the same volume moved (127 m 3 = 10 passes for forwarding and 45 for skidding), due to the different workload capacities of the two extraction methods being examined.

Physical soil parameters: BD and PR
Our findings confirmed the general rule that the passage of a loaded ground-based machine causes soil compaction (in terms of both BD and PR increase) along the trails (Williamson and Neilsen 2003;Han et al. 2009;Cambi et al. 2015).
Indeed, BD increased for both forwarding and skidding in InR after 10 passes, driving in both uphill and downhill directions. This difference is attributed to the direct (under tyre in InR) effect of the pressure exerted on the soil. Nonetheless, a significant increase of BD was found also between undisturbed area and in BtwR, because of an indirect (between left and right tyres in BtwR) effect of the pressure exerted on the soil. These findings confirm that the impacts of machine traffic on soil may not be limited only to the contact area of tyres on terrain, but also on the surroundings (Solgi et al. 2016). In fact, Solgi et al. (2016) demonstrated that the soil compaction expressed as BD occurred as far as 1 m from the wheel ruts and increased with slope gradient and traffic intensity, which explains the differences in our results between undisturbed area and BtwR in terms of BD.
In forwarding we found a significant difference between InR and BtwR in terms of BD, with the highest impacts in InR. However, the difference between InR and BtwR is not noticeable in skidding. This distinction between skidding and forwarding may be attributable to the different forces involved in the wheel-soil interaction process and the differences in load distribution for each extraction method. During forwarding, indeed, the load is carried on the trailer and the impact on soil is caused by the pressure exerted by the machine tyres combined with wheel slippage on the ground. During skidding, in contrast, the bottom ends of the logs are lifted and connected to the winch and the top ends are dragged along the ground-thus, the impact on soil is caused by a combination of the pressure exerted by the tyres and the pressure/scarifying effect of the top ends of the logs. When the tyre pressure exceeds the soil bearing capacity-especially under increasing tractive demand-subsequent wheel slippage can induce pronounced shearing processes at the soil surface (Edlund et al. 2013), thereby increasing soil compaction (Eliasson and Wästerlund 2007). Therefore, although this is the case for both skidding and forwarding, during skidding the logs dragged along the skid trail exerted pressure and scarified the soil surface, thus affecting compaction, soil mixing and displacement. The significant difference between InR and BtwR found for forwarding was not recorded in skidding due to the pressure exerted by the top-ends of dragged logs in contact with the ground. Accordingly, in BtwR, it was expected that higher BDs and PRs would be found for skidding than for forwarding, due to the effect of the dragged logs on the ground. However, 10 passes of dragged logs were not sufficient to determine a significantly higher compaction in comparison with forwarding, in either downhill or uphill directions. On the contrary, after 127 m 3 of wood volume moved the soil compaction was affected by the driving directions (Table 4). Nevertheless, results refer to the specific work conditions applied in this study; different lengths of logs and different distances of bottom ends from the ground could change load distribution varying the forces applied on the ground.
It is interesting to highlight that, although BD values were always higher after machines pass, PR showed no significant differences in BtwR for forwarding-both downhill and uphill-as well as for skidding downhill. However, PR increase for skidding downhill in BtwR became significant increasing the number of passes, i.e. considering the same wood volume moved as that of forwarding (127 m 3 ). Considering these results, it appears that PR in BtwR is less affected than BD by the indirect effects of tyre pressure. In fact, the analysis of variance showed a strong influence of the wood extraction methods on BD (Table 4), difference that decreases with the increasing number of machine passes (10 forwarding and 45 skidding passes to obtain 127 m 3 wood volume moved).
After 127 m 3 transported, similar values of BD between skidding and forwarding were registered in uphill extraction, while in downhill operation the effects of skidding were opposite to those of forwarding considering InR/BtwR; regarding PR, the trend seems similar, but differences were not as strong for BD.
Differences in BD levels of between uphill and downhill skidding may be explained by an uneven load distribution on the machine's axles, depending on slope direction (Jamshidi et al. 2008;Jourgholami et al. 2014) and by the increased wheel slippage and vibration encountered when skidding uphill compared to downhill. In this context, in the downhill direction, the machine's weight is more homogeneously distributed between front and rear axles. Furthermore, the weight of the logs sustained by the winch does not reduce the load on the front axle as occurs in the uphill direction. Therefore, the different distribution of weights between uphill and downhill directions implies a different application of traction force by the tractor on the ground, varying both compaction and shear effects on soils due to the threedimensional transmission of forces (Horn et al. 2007). For these reasons, in some cases the effect of dragged logs on the ground may be greater than tyre pressure, thus explaining higher BD values in BtwR than in InR in the downhill direction after skidding 127 m 3 . On the contrary, in the uphill direction, the tractor weight distribution is mainly on the rear axle, and the weight of logs acting on the winch increases the rear-up effect, thus concentrating the pressure of the rear axle tyres on the ground (McCallum B 1993). Moreover, the greater the slippage the greater the exposure of mineral subsoil, which naturally has a higher BD than the surface layer (Jourgholami et al. 2014), further underlining the reasons why greater BD resulted in the uphill direction, as also reported by Sidle and Drlica (Sidler and Drlica 1981). Finally, the examined machines-tractor with winch or trailer-may not be the optimal solution available in comparison with machines specifically built for forest operations (i.e. forwarder and skidder). In fact, the ground tyre pressure for the loaded trailer, 11,772 kg in total, was high (173 kPa on average per tyre) due to the narrow and small tyres. This could negatively affect the severity of damages on soil, that may be reduced using a skidder or a forwarder, normally equipped with larger tyres with a low pressure of inflation.

Rut measurements
The analysis of rutting after machine traffic was focused on two main aspects: (i) identifying the best method for rutting data collection and processing and (ii) assessing the effects of machine passes in terms of soil volume moved.

Comparison of DSMs obtained by SfM vs. PLS
The other aspect investigated in this study was rutting formation, analysing in depth the effects of tractors with winch or trailer in terms of both rut (volume reduced: VR) and bulge (increased volume: VI). The two methods applied and compared-PLS and SfM-have not shown significant differences in terms of results. The comparison of DSMs obtained by PLS measurements and SfM methods showed small differences (i.e. 97% of the plot area with difference < 2 cm). However, the use of more precise equipment, such as a total station, can further increase the precision of soil reconstruction in both methodologies (Pierzchała et al. 2016;Starke et al. 2020). SfM is the preferred method for technical reasons; indeed, in comparison with PLS, photogrammetry has much higher resolution (almost double the point m −2 in the dense cloud) and spatial distribution of the point cloud. In fact, the point resolution of laser scanning clouds increases for areas closer to the PLS, whereas it decreases at longer distances, which is one of the main limitations of this technology (Nadal-Romero et al. 2015). It is therefore recommended to minimize the scanning distance when possible (Heritage and Hetherington 2007). In contrast, the resolution of multiple views with the same imaging height and overlap used in SfM were not affected by object distance (Smith and Vericat 2015). Moreover, when using PLS, the data check operation must be performed afterward (i.e. not immediately in the field). In photogrammetry, however, it is possible to check on the camera screen if the collected images are of good quality and if they correctly include the GCPs, covering the studied plot as planned, thus easily and promptly identifying and correcting errors during data collection, avoiding further field surveys in case of errors as may happen when using PLS. However, two of the main problems related to SfM are object occlusion and reflection. In the forest, vegetation (grass and/or shrubs) can cover the ground and water can reflect the signal in the rut, both compromising the accuracy of the method (Pierzchała et al. 2016), implying the necessity of choosing almost dry conditions for data collection. At the same time, one of the main problems in modelling a terrain based on a point cloud obtained by PLS is distinguishing points reflected by the terrain from those reflected by other objects (Koreň et al. 2015). In our work, we used a simple and effective filtration algorithm that removed all points located at a certain height above the ground (Koreň et al. 2015). However, the DSMs after machine traffic could be also affected by the presence of harvesting residues on the forest floor or stones overturned by skidding. In this case, a visual detection of such objects is more difficult in PLS point clouds in comparison with coloured photogrammetry point clouds' negative effects on the reliability and accuracy of the method. Nonetheless, logging residues left on the ground should be removed from the sampling area before data collection.
In conclusion, both PLS and SfM photogrammetry have increased the capacity to provide large data sets for the estimation of terrain modifications; without these technologies, the common method was manual, which was time-consuming and less reliable than PLS and SfM (Pierzchała et al. 2016;Haas et al. 2016;Marra et al. 2018).

Quantification of rutting
Photogrammetry showed interesting results on the effects of skidding and forwarding as causes of rutting. In fact, results showed that after 127 m 3 of wood volume had been moved, the tractor with the winch caused ruts similar to those caused by the tractor with the trailer in terms of soil volume reduced (VR). On the contrary, considering the increased volume (VI: bulges), skidding caused an effect four times greater than that of forwarding, with 7.36 m 3 against 1.68 m 3 of VI, respectively, per 100 m of trail. The reason for this considerable difference is the effect of log top ends dragged on the ground, and it is partially explained by what has already been discussed regarding the effects on soil physical parameters (see Discussion -physical soil parameters: BD and PR). In fact, rutting is caused by machine wheels, both in forwarding and skidding-but in skidding, the passing of the end of the dragged logs may change the impacts on soil (Wood et al. 2003;Cambi et al. 2018a). In detail, the ends of the dragged logs scratch and displace a certain quantity of soil in the dragging direction during each extraction trip. In this way, the soil is reshaped after the passage of tractor and the ruts left by the wheels are hidden. This happens especially in slope changes when the logs move and rearrange (Heninger et al. 2002;Williamson and Neilsen 2003;Agherkakli et al. 2010;Cambi et al. 2018a). Moreover, in skidding wheels cause a higher soil displacement than in forwarding due to pulling action. Looking at DSMs, rutting values are visually confirmed; in fact, in forwarding, ruts were clearly identifiable, whereas bulges along the trail were not visible. In contrast, ruts were not clearly evident on skidding trails, while bulges were. This difference negatively affected skidding in terms of risk of soil erosion, being the volume of soil moved higher than forwarding. In fact, the soil displaced by dragged logs and disturbed due to the traffic can be very vulnerable to erosion (Bagheri et al. 2013), causing negative indirect impacts on forest ecosystem, water quality and land stability (Frey et al. 2009); in this way, the measured volumes of bulges can be considered as a reference value of the potential entity of short-time erosion along the trail.

Conclusions
The main aim of this study was to describe the effects of wood extraction on both the physical properties of soil and the rutting caused by tractors with winch or trailer in different wood extraction directions. Our findings highlighted that the use of two different wood extraction methods can affect the soil impact differently in opposing driving directions. Within the limits of our experimental conditions, the study results highlighted that: i. the role of the driving direction was found to be irrelevant when operating a four-wheel driven tractor and a trailer with one of its two axles driven. On the contrary, in order to reduce soil compaction, downhill skidding should be preferred to uphill skidding; ii. comparing skidding and forwarding, downhill skidding caused lower soil compaction than forwarding, while no differences were found uphill. For this reason, in downhill extraction skidding should be preferred to forwarding, while uphill-being soil compaction comparable-the choice can be made on the basis of other parameters, such as technical and logistics issues; iii. skidding has a higher effect on soil displacement than forwarding; In addition, our study applied and compared the performance of portable laser scanning and image-based models derived by SfM photogrammetry to evaluate soil rutting. Indeed, the use of both portable laser scanning and imagebased models derived via SfM photogrammetry were accurate methods for the creation of high-resolution DEMs and may be very useful to assess the impact of forest operations on soil. However, considering the logistics of data collection in forestry, photogrammetry is likely the best solution.
Finally, the results of this study are an important piece of knowledge aiming to reduce negative effects on soil due to ground-based extraction methods in forest operations. Following the principles of Sustainable Forest Operations (Marchi et al. 2018), modern logging operations must minimize these impacts, and this results can be useful to plan effective mitigation strategies reducing negative effects of logging on soil.
Funding This research was supported by SKIDDFORW project, funded by the University of Florence. Open access funding provided by Università degli Studi di Palermo within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.