Comparison of Selected Terramechanical Test Procedures and Cartographic Indices to Predict Rutting Caused by Machine Traffic during a Cut-to-Length Thinning Operation

Abstract

Timber harvesting operations using heavy forest machinery frequently results in severe soil compaction and displacement, threatening sustainable forest management. An accurate prediction of trafficability, considering actual operating conditions, minimizes these impacts and can be facilitated by various predictive tools. Within this study, we validated the accuracy of four terramechanical parameters, including Cone Index (MPa, Penetrologger), penetration depth (cm, Penetrologger), cone penetration (cm blow⁻¹, dual-mass dynamic cone penetrometer) and shear strength (kPa, vane meter), and additionally two cartographic indices (topographic wetness index and depth-to-water). Measurements applying the four terramechanical approaches were performed at 47 transects along newly assigned machine operating trails in two broadleaved dominated mixed stands. After the CTL thinning operation was completed, measurement results and cartographic indices were correlated against rut depth. Under the rather dry soil conditions (29 ± 9 vol%), total rut depth ranged between 2.2 and 11.6 cm, and was clearly predicted by rut depth after a single pass of the harvester, which was used for further validations. The results indicated the easy-to-measure penetration depth as the most accurate approach to predict rut depth, considering coefficients of correlation (r_P = 0.44). Moreover, cone penetration (r_P = 0.34) provided reliable results. Surprisingly, no response between rut depth and Cone Index was observed, although it is commonly used to assess trafficability. The relatively low moisture conditions probably inhibited a correlation between rutting and moisture content. Consistently, cartographic indices could not be used to predict rutting. Rut depth after the harvester pass was a reliable predictor for total rut depth after 2–5 passes (r_P = 0.50). Rarely used parameters, such as cone penetration or shear strength, outcompeted the highly reputed Cone Index, emphasizing further investigations of applied tools.

1. Introduction

Currently, the forest in Germany and other countries suffer from several climate change-induced stressors like drought, fires and other extreme weather conditions [1,2]. One consequence of these weather phenomena is an ongoing large scale bark beetle infestation [3]. Thus, prompt salvaging of infected trees is necessary for forest protection [4,5,6] and to limit economic losses for forest owners. In some regions, the required machine-operating trails for off-road traffic still have to be established to ensure the necessary year-round access to the affected stands. Despite this, off-road traffic is currently a general matter of debate, since it causes multiple adverse impacts to soil characteristics [7,8,9,10,11], and thus the stand growth (e.g., [12]). Negative impacts like compaction (e.g., [13]) and soil displacement (e.g., [14]) are the consequence of heavy machine traffic on forest soils. Following Horn et al. [15], soil displacement, caused by shearing forces and soil compression under wet conditions, leads to ruts. Soil moisture content, terrain slope, soil properties and the operating vehicle are decisive for the creation of ruts (e.g., [16]). Additionally, the cumulative weight associated with the number of passes is linked to rut formation [7,17,18]. Since ruts show a strong relationship with decreased forest productivity in the affected areas [19], it is a priority of sustainable operation management to keep their formation to a minimum.

Estimating soil trafficability is a key aspect for preserving soil integrity and thus the forest as an ecosystem on the one hand, and for optimizing timber harvest and efficient machine utilization on the other. In theory, soil trafficability can be derived from terrain characteristics, machine configurations, and weather conditions. An accurate prediction of rutting is necessary for sustainable, low-impact harvesting operations, since it supports the mitigation and avoidance of impacts [18,20,21] on machine operating trails, which are established without a pavement or any other fixation. The assessment of forest sites’ trafficability continues to be commonly based on the appraisal of the forester in charge (e.g., [22]), mainly requiring a long-year experience to enable a sufficient estimation of current harvesting conditions in a designated area. However, because of varying weather conditions considerably affecting the trafficability, the assessment of current conditions and subsequently the timing of operations is a challenging task. In many cases, the machine operator himself is faced with this task when arriving at the site of operation. Usually, the machine operator starts the operation, until the occurring ruts are individually judged as too severe. Consequently, the operation has to be stopped, requiring a cost-intensive relocation of the machine to alternative stands with better trafficability. Hence, lacking knowledge of current trafficability on stand levels may result in manifold drawbacks such as financial shortcuts and, even worse, ecological damages with long lasting, severe impacts on tree growth, water infiltration, soil organisms and soil erosion [23].

A prediction of such damages can be facilitated by terramechanical test procedures (e.g., [24]) or by trafficability maps [25]. Recent approaches show the possibility to incorporate information on trafficability, measured during the first machine pass, in such maps [18,26,27]. Both terramechanical tests and trafficability maps can sufficiently support the operational flexibility of company management and entrepreneurs as well as machine operators, in order to avoid and reduce disturbances to forest soils, with associated high costs. Some of the equipment used for assessing soil’s bearing capacity has already been validated towards the predictability of ruts. For example, Sirén et al. [18,28] reported a correlation between rut depth after harvest and Cone Index (CI). The same index has been frequently used for validations of accessibility maps, where CI was used instead of assessing actual machine traffic (e.g., [29]). Less commonly used, but already reported by Heubaum et al. [30], a vane meter can be used to assess soil trafficability through the determination of shear strength. Especially if measured in a mineral soil depth between 10 and 15 cm, this rather simple instrument has been reported as promising regarding the accuracy of prediction. Besides, the relatively unknown dual-mass dynamic cone penetrometer has been used in the estimation of pavements moduli [31,32]. However, as far as we know, it was not validated in forest soils up to now. Certainly, numerous studies revealed the link between soil moisture content and rutting [23,33], giving this parameter a high relevance for the prediction of rutting.

Subsequently, two pivotal approaches of remote analyses for the simulation of soil moisture content were chosen in addition. The calculated indices topographic wetness index (TWI) [34,35] and depth-to-water (DTW) index [36,37,38] link the terrain characteristics with the soil wetness by the simulation of accumulating water runoff. In principal, both indices quantify the contributing area to a distinct cell within a grid, and calculate the likelihood for cells to be water-saturated or moist. Especially in dense stands, with a lacking network of machine operating trails, visibility is restricted for machine operators. The cartographic knowledge of moist and therefore sensitive areas, in particular when available through geo-referenced digital maps on the forest machine’s on-board computer screen, can highly contribute to the avoidance of deep ruts and other soil impacts [39].

Still, the above-mentioned approaches for the prediction of rutting differ in regard to accuracy, financial efforts, handling and applicability, leading to specific advantages and drawbacks. Nevertheless, comparative studies between different approaches are lacking. In order to fill this knowledge gap, we conducted a study which aimed to assess relevant terramechanical test procedures and cartographic indices, all assumed to predict rutting during forest operations. In particular, the following questions were addressed:

• Can rut depth after a cut-to-length (CTL) thinning-operation be predicted by the two used cartographic indices (i.e., depth-to-water and topographic wetness index)?
• Which of the applied terramechanical test procedures show a reliable response with occurring rut depth, in terms of a high Pearson coefficient of correlation?
• Is rut depth formation after the first machine pass (i.e., facilitated by a harvester) a reliable figure to predict total rutting of a consecutive CTL thinning operation?

To answer these questions, a field trial was conducted at newly assigned machine operating trails. There, cartographic indices were calculated and in-situ measurements were performed in advance to a regularly planned CTL thinning operation. Rut depth was estimated after the first and final machine pass and analyzed in relation to cartographic indices and previously conducted field measurements.

2. Materials and Methods

2.1. Stand Characteristics

The study was performed in two stands, both located near the town of Uslar in Lower Saxony, Germany. The larger stand B (Figure 1B, area: 7.7 ha, x: 9.61, y: 51.71, WGS 84) was characterized by 95-year-old beech (Fagus sylvatica), showing a mean diameter at breast height of 28 cm. This stand was intermixed with grouped Norway Spruce (Picea abies, age: 75 years). The smaller stand A (Figure 1A, area: 1.9 ha, x: 9.60, y: 51.70, WGS 84) was stocked by 40-year-old beech, scattered with European larch (Larix decidua) and Norway Spruce. The stands were growing on Cambisol [40], based on a silicate bedrock. The terrain was smooth with a slight slope, less than 30%. The long-term mean of annual precipitation amounts 892 mm [41]. Yet, overall dry conditions were present during the field campaigns, caused by a low precipitation of 345 mm between March and August [42].

2.2. Harvesting Operation

A CTL thinning operation was performed at the study sites. Felling and processing were executed by a single grip harvester (see

Table 1. Machine specifications of the used 8-wheel machines, a harvester (Ponsse Bear) and a forwarder (Ponsse Buffalo), conducting the investigated CTL thinning operation.

Character	Unit	Ponsse Bear	Ponsse Buffalo
power	kW	260	210
typical mass	Mg	24.5	19.8
loading capacity	Mg	-	15.0
tire type		Nokian Forest King TRS 2	Alliance Forestar 344, 20 PR
tire size (width, diameter)	mm	750, 1485	710, 1340
inflation pressure	kPa	600	500

for machine specifications) on 4 August 2020 (

Table 2. Overview and date of performed in-situ measurements. For the study, volumetric and gravimetric soil moisture content (SMC_VOL, vol%; SMC_GRAV, % (M_FRESH and M_DRY, Equation (5)), respectively), commonly used and a modified Cone Index (CI, CI_MOD, respectively, MPa), penetration depth (PD, cm), penetration resistance by means of dual-mass dynamic cone penetrometer (DCP, cm blow⁻¹) shear strength (τ, kPa) initial and post-operational soil bulk density (SBD_INIT, SBD_POST, respectively, g cm⁻³) were measured. For the estimation of rut depth increments, transect profiles were measured (Dz_init, Dz_H, Dz_F, Figure 3, cm).

Date (2020)		Objective	Measurement
15 July	initial measurements of 90 transects	reference profile to estimate rut depth increment and	Dzinit
		terramechanical parameters	SMCVOL, CI, CIMOD, PD, DCP, τ, SBDINIT, SMCGRAV
4 August	harvester performed the felling and processing
10 August	measurement of profiles	rut depth after harvester	DzH
17–21 August	forwarder excavated timber
1 September	post-operational measurements	rut depth after forwarder, total rut depth and	DzF
	on the remaining 47 transects	post-operational soil bulk density and moisture content	SBDPOST, SMCGRAV

). In the days between felling and forwarding (Table 1, 17–21 August), 45 mm of precipitation occurred at the site. On average, 36 m³ ha⁻¹ were harvested at the study area by clearing the machine-operating trails and thinning the forest stands.

The operating harvester resulted in a single machine pass on the machine operating trails (number of machine pass = 1). The subsequently operating forwarder was equipped with a GNSS device, constantly tracking the machine position during timber extraction. Thus, the number of passes, as illustrated in Figure 1, was assigned to all 47 measuring transects.

2.3. Field Measurements

2.3.1. Rutting on Machine-Operating Trails

Following the concept of permanent machine operating trails, a new system of parallel running machine-operating trails, with a spacing of 24 m from the center line, was established in both stands to make the area of 9.6 ha in total accessible for forest machines. Therefore, the projected machine operating trails were marked in the stands. On a total of 19 newly planned machine operating trails, 90 perpendicularly positioned transects with a length of 4 m each, were created. At both ends, wooden pegs were placed, in order to mark each transect and to position a reference beam (Figure 2), used to measure the rutting process during the conducted harvesting operation. For that purpose, the initial profiles along each transect were measured using a yardstick at 20 cm spacing.

This initial profile or reference surface was defined by the distance between the beam (levelled at the wooden pegs) and the surface, giving Dz_init (Figure 3). This (1.) measurement was done after loose humus material had been removed. Profiles along each transect were repeatedly measured: (2.) after the harvester passed, giving Dz_H (Figure 3) and (3.) after the timber extraction via forwarder was accomplished, giving Dz_F (Figure 3). Several wooden pegs, used to position the reference beam, were severely displaced during the harvesting operation, reducing the initial 90 measuring transects to the remaining 47 ones (Figure 1), since the rut depth could not be quantified anymore. When Dz_H and Dz_F were measured, the position of the tracks on the transect was captured as a dummy variable (t, where a visually detectable track was defined as “L” (left machine track) or “R” (right machine track,

Table A1. Example of a measured transect, used to calculate rut_H, rut_F and rut_T (Section 2.3.1), by means of the code shown in Appendix B. Respectively to Dz_H and Dz_F, t_H and t_F were noted, to capture the position of visible machine tracks across the measuring transect.

transect	tr.length	DzINIT	DzH	DzF	tH	tF
41	−200	36	34	34	0	0
41	−180	32	33	33	0	0
41	−160	32	32	32	0	0
41	−140	32	32	31	0	0
41	−120	30	32	31	L	0
41	−100	31	33	33	L	L
41	−80	31	32	33	L	L
41	−60	29	30	35	L	L
41	−40	30	29	32	0	L
41	−20	27	28	29	0	L
41	0	29	29	28	0	0
41	20	29	30	29	0	0
41	40	30	29	30	0	0
41	60	29	30	29	0	0
41	80	31	32	31	0	0
41	100	29	32	31	R	0
41	120	32	33	34	R	R
41	140	33	34	36	R	R
41	160	35	36	37	R	R
41	180	37	37	37	0	R
41	200	38	39	38	0	R

), used for the calculation of occurring rutting (Appendix B). Consequently, where t equals “L” or “R”, the maximum difference between Dz_H and Dz_init (Figure 3) of the left and right track were used to calculate the rut depth after the operating harvester, rut_H (cm), according to Equation (1):

$${\mathrm{rut}}_{\mathrm{H}}={\displaystyle \sum }_{}^{n=2} \frac{D{z}_{H,i}-D{z}_{init,H,i}}{n}$$

(1)

Used variables for Equations (1)–(3) are illustrated in Figure 3, where i represents “L” or “R”; the applied R-grammar is shown in Appendix B.

Following the same approach, rutting after forwarding (rut_F, cm) was calculated according to Equation (2), whereas the differences between Dz_F and Dz_H were used to quantify rut_F, describing the additional rutting after the harvester has passed already:

$${\mathrm{rut}}_{\mathrm{F}}={\displaystyle \sum }_{}^{n=2} \frac{D{z}_{F,i}-D{z}_{H,F,i}}{n}$$

(2)

Since total rut depth (rut_T, cm) cannot be calculated as the sum of rut_H and rut_F, as both machines might not drive exactly in the same machine track, the same procedure was applied to quantify rut_T. The differences between either Dz_F or Dz_H, and Dz_init were analyzed on tracks, after the CTL thinning operation was accomplished, as illustrated in Figure 3, according to Equation (3):

$${\mathrm{rut}}_{\mathrm{T}}={\displaystyle \sum }_{}^{n=2} \frac{D{z}_{F|H,i}-D{z}_{init,F,i}}{n}$$

(3)

2.3.2. Soil Samples

Along the newly assigned and still uncleared machine operating trails, soil samples were collected at randomly chosen measuring transects. Pre- and post-operational sampling was done at 15 measuring transects (Figure 1, rhombs), out of the 90 initial transects. Samples collected at nine of these 15 transects could be used to correlate parameters derived from soil samples with occurring rut depth, due to the above-mentioned reduction of measuring transects to 47 remaining ones (Section 2.3.1). Still, the six remaining sampling positions (Figure 1, black rhombs) were used for correlations with terramechanical test procedures only (Section 2.3.3).

The soil samples were either collected at the potential right or left machine track (position was assumed during the initial sampling), using 100 cm³ sampling rings. Pre- and post-operational samples were taken in a depth between 10 and 15 cm of the mineral soil. Subsequently, initial soil bulk density (SBD_INIT, g cm⁻³), as well as post-operational soil bulk density (SBD_POST, g cm⁻³) were used to calculate compaction (comp, %), according to Equation (4):

$$comp= \frac{SB{D}_{POST}- SB{D}_{INIT}}{SB{D}_{INIT}}\ast 100$$

(4)

The fresh samples were weighed, giving an initial mass (M_FRESH, g), then dried in an oven (105 °C) until mass constancy was reached. Afterwards, samples were weighed again, giving M_DRY (g), used to calculate gravimetric soil moisture content (SMC_GRAV, %), according to Equation (5). Additionally, soil texture was analyzed on three mixture probes, containing all soil samples, according to the approach of Durner et al. [43]

$$SW{C}_{GRAV}= \frac{M_{FRESH}- M_{DRY}}{M_{DRY}}\ast 100$$

(5)

2.3.3. Terramechanical Test Procedures

Despite SBD_J and M_J, both being measured before and after the CTL thinning operation, the remaining soil parameters were measured in advance of the operation only (Table 2), describing initial conditions. These initial conditions were further used to address the posed research question concerning the comparison among the included terramechanical test procedures. For that, initial measurement results were associated with rut_H. Soil moisture content, penetration resistance, depth and shear strength were measured on the 90 initial measuring transects, but only 47 of them could be used to be correlated with rut_H. The measurements by the used instruments (Figure 4) were spaced by 10 cm to each other, with a spacing of 50 cm along each transect (Figure 2):

• Moisture meter: Volumetric soil moisture content (SMC_VOL, vol%) was quantified in the mineral topsoil, where a 57 mm long TDR probe (HH2-moisture meter, Delta-T-Devices, Cambridge, UK) was inserted from above, after the removal of humus. This moisture meter measures volumetric moisture content, θv, by responding to changes in the apparent dielectric constant of moist soil, resulting in a ratio between the volume of water and the total volume of the soil sample [44]. Seven measurements on each transect were averaged, giving SMC_VOL.
• Penetrologger: Penetration resistance was measured, using a handheld Penetrologger (1.0 cm², 60° cone, Eijkelkamp Soil and Water, Giesbeek, The Netherlands). This device captures the soil penetration resistance for each centimeter by means of a load cell, whereas values of the upmost 15 cm of mineral soil were averaged giving a Cone Index (CI, MPa), as mean value for seven penetrations on each transect. Based on previous experience by the authors [45], a modified Cone Index (CI_MOD, MPa) was calculated in addition. In contrast to the estimation of CI, soil penetration values between 10 and 20 cm were considered to quantify CI_MOD, due to the high variance of penetration resistance in the upmost centimeters. Besides, the total penetration depth, captured by the Penetrologger, was averaged for each measuring transect, giving PD (cm).
• Dual-mass dynamic cone penetrometer: Since we decided to keep the time demand for measurements approximately similar between the compared methods, the number of samples was reduced for this instrument. Consequently, one measurement, consisting of six hammer blows, was done in the middle of each transect. The incremental penetration depth was captured and used to calculate the corresponding parameter, derived by the dual-mass dynamic cone penetrometer, DCP (cm blow⁻¹), defined as the average of penetration depth per hammer blow, until it exceeded 15 cm penetration depth.
• Vane tester: The used Eijkelkamp (Eijkelkamp Soil and Water, Giesbeek, The Netherlands) field inspection vane tester is an instrument for in-situ measurements of shear strength through vanes of different sizes. Measurements are conducted through a spiral spring, detecting the torque, which needs to be applied to a handle in order to displace the soil through the vane. The used (and smallest, 16 mm × 32 mm sized) vane allows to cover readings up to 260 kPa, by an accuracy within 10% [46]. The current measurements were done in a mineral soil depth of 10 to 15 cm, as recommended by Heubaum et al. [30], giving shear strength (τ, kPa) as mean value for each measuring transect.

2.4. Trafficability Maps

During the field campaigns, multiple positioning measurements were conducted using a handheld GPS device (Oregon 700, Garmin Ltd., Olathe, KS, USA). Later, these waypoints were averaged to get the nearest approximation of the transects’ position. From the 90 initial measuring transects, 47 were used for the correlation between measured rut depth and calculated indices. These point positions were used to join field measured data with remote processed data using QGIS (version 3.10.10, Open Source Geospatial Foundation Project [47]). Two trafficability indices for the identification of wet and therefore sensitive areas within forest stands were calculated using QGIS.

2.4.1. Depth-to-Water Index

Digital terrain models (DTM) are the only prerequisite input data to calculate depth-to-water (DTW) maps. Within this study, a DTM with a grid size of 1 m has been used, originally collected by the State Office for Geoinformation of Lower Saxony [48]. Based on this model, DTW maps were created according to Murphy et al. [36], as extensively described already [25,36,37,49,50,51,52,53]. In short, a flow accumulation model of the DTM is generated using the D8 flow algorithm. When the flow accumulation of a cell reaches a certain value, the very cell is defined as the starting point of a flow line, which continues the flow downhill until leaving the mapped area. The required area to start a flow line is defined as “Flow Initiation Area” and herein was set to 4 ha, as suggested by Jones and Arp [29] for likewise conditions. Subsequently, a least-cost function calculates the least slope path from each cell within the grid towards the nearest flow line [36,37]. Hence, low values of the metric DTW index indicates a spatial proximity (in a vertical direction) to the nearest flow line, with assumed wet or poorly to imperfectly drained areas. A DTW index of up to 1 m is usually assumed to show areas with water-saturation, which are associated with a high susceptibility towards rutting. Contrarily, high values of the index generally indicate dryer and therefore accessible areas.

2.4.2. Topographic Wetness Index

In addition, the topographic wetness index (TWI), commonly used to quantify topographic control on hydrological processes was calculated according to Sörensen et al. [54]. This index combines the local upslope contributing area (a) and slope (β), given by Equation (6):

$$TWI=\ln \frac{a}{\tan \beta }$$

(6)

Within this study, we aggregated the grid size of 1 m to receive a grid of 5 m resolution [53]. The procedure Zevenberg Thorne was chosen within a GDAL algorithm to calculate the slope grid, since the terrain was smooth. To generate the required flow accumulation grid, a SAGA algorithm was applied, which filled occurring sinks during reprocessing. Subsequently, having slope and contributing area given by the flow accumulation, TWI could be calculated according to Equation (6).

2.5. Data Analysis

All data were merged and analyzed using R core (version 4.0.2, R Foundation for Statistical Computing, Vienna, Austria [55]), interfaced with R studio (version 1.3.959, PBC, Boston, USA [56]). Normal distribution of CI, CI_MOD, DCP, τ, rut_H and rut_F was assessed and approved via QQ-plots. Both cartographic indices were adjusted, DTW was transformed into a dummy-variable, with a margin of 1 m, TWI was log-transformed. Individual linear models were fitted between parameters. Coherently, Pearson coefficients of correlation (r_P) were used to evaluate the quality of the given models, and were depicted applying a modified correlation matrix [57]. Residuals of linear models were visually assessed to identify outliers. In addition, values from transects with exceptionally high traffic were identified. Accordingly, four observations were removed. Afterwards, normal distribution of residuals was tested and confirmed by Shapiro–Wilk tests. Initial gravimetric soil moisture content was compared with post-operational, by means of a one-sided paired t-test. An unpaired t-test was applied to test for differences of rut_H measured within predicted sensitive areas (DTW) and remaining ones. Analyses were usually performed on 47 observations (90 initial transects, 43 of them unusable, since rut depth could not be measured, or seemed unreliable according to the outlier analysis), or 15 observations for variables derived from soil samples. The correlation between soil bulk density and rut_H was performed on nine remaining observations. All tests were done using a significance level of 0.05, tendencies were reported up to a significance level of 0.10. Values within the text are usually given as mean ± standard deviation.

3. Results

3.1. Soil Properties

The Cambisol [40] at the study area, composed of 33 ± 5% sand, 35 ± 10% silt and 32 ± 7% clay, showed a mean penetration depth of 31 ± 15 cm, derived from measurements using the Penetrologger. Across the study site, the initial volumetric moisture content averaged at 29 ± 9 vol% (spatial distribution is depicted in Figure 5), before the operation had been started. Initial gravimetric soil moisture content (SMC_GRAV), estimated on 15 gathered soil samples, described similar values, with an average of 31 ± 6%. Due to 45 mm of rainfall, occurring between the time of felling (harvester) and forwarding (Table 1 and Table 2), SMC_GRAV slightly increased to 34 ± 7% (p = 0.046). Apparently, the increase of SMC_GRAV was higher on measuring transects with low initial SMC_GRAV (Figure 6A), compared to rather moist initial conditions. Soil bulk density across the study site was quantified with 1.2 ± 0.1 g cm⁻³. Through the harvesting operation, this initial density increased by 10% (p = 0.017), up to 1.3 ± 0.1 g cm⁻³. It was observed that measured compaction showed higher extents on measuring transects with lower initial soil bulk density (Figure 6B).

3.2. Rutting

The clearing of the 19 machine operating trails and the thinning of the adjacent stand area (Figure 1) was conducted by the 24.5 Mg heavy harvester (Table 1), resulting in on average 3.6 ± 1.7 cm deep ruts (rut_H), ranging from −0.4 to 7.0 cm. The subsequently driving forwarder, trafficking the same machine operating trails, resulted in additional average rut_F of 3.8 ± 2.1 cm, ranging from 0.0 to 9.9 cm. Although differences between rut_H and rut_F were low, rut_F clearly explained more variance of total rut depth (rut_T), with r² of 0.63, compared to r² of 0.25, when rut_H was correlated with rut_T. Still, rut_T responded to rut_H, given as rut depth after the first machine pass (p < 0.001), as can be seen in Figure 7A. Thereby, rut_T could be quantified with 6.3 ± 2.1 cm (2.2 to 11.6 cm). Under the given operational conditions of the current study, rut_T showed no dependency of traffic frequency (2–5 passes, Figure 7B).

3.3. Correlations with Terramechanical Tests

Terramechanical tests were applied on every measuring transect along the previously marked machine operating trails, to capture initial conditions. Descriptive statistics for these tests are given in

Table 3. Descriptive statistics of terramechanical parameters, measured prior to harvesting operation. Parameters are regular and modified Cone Index (CI, CI_MOD, respectively) and penetration depth (PD), measured by Penetrologger, DCP (cm blow⁻¹, dual-mass cone penetrometer) and shear strength (τ, kPa, vane tester). With the number of observations, mean value, standard deviation, extreme values, first and third quantiles.

Parameter	n	Mean	SD	Min.	0.25	0.75	Max.
CI (MPa)	47	1.52	0.31	1.01	1.29	1.70	2.53
CIMOD (MPa)	47	2.24	0.61	1.25	1.78	2.51	3.61
PD (cm)	47	30.86	14.78	10.29	18.00	41.00	62.86
DCP (cm blow−1)	47	4.18	4.22	0.00	1.77	4.94	21.60
τ (kPa)	47	214.02	69.63	127.43	150.14	265.57	365.14

.

Due to high uncertainties regarding the load mass of the forwarder, we decided to use rut_H to validate the surveyed cartographic indices, and terramechanical test procedures. Three of the contained parameters, CI, CI_MOD and shear strength (τ), are capturing forces, required to overcome a given resistance within the soil. Analogously, these parameters showed inverse responses to rut_H, as can be seen in Figure 8 and Section 4.2. Whereas a correlation can be assumed for τ (p = 0.088), an unreliable response was observed for CI (p = 0.415). Regardless of a standard definition of CI, given as the mean penetration resistance of the upmost 15 cm soil, we also generated a modified CI. For this, penetration resistance of a soil depth between 10 and 20 cm was averaged, giving CI_MOD, which was significantly linked with rut_H (r_P = −0.30, p = 0.040, Section 4.2). In addition to CI and CI_MOD, penetration depth (PD) was also derived from measurements by the Penetrologger. Whereas CI possessed an insufficient prediction of rut_H, the latter showed a significant correlation to PD (p = 0.001) with a high r_P of 0.44. As indicated in Figure 8, DCP, measured by means of the dual-mass cone penetrometer, correlated with rut_H (p = 0.020), resulting in the second highest value of r_P (0.34).

3.4. Correlations with Cartographic Indices

A DTW map with a flow initiation area of 4 ha was calculated for the study area. There, the values of DTW ranged between 0 and 8.7 m. The overall mean at the analyzed transects was 3.6 ± 2.3 m. Out of those, four measuring transects exhibited a DTW index below 1 m (

Table 4. Mean values and standard deviations (SD) of rut depth after the harvester (rut_H, Table 1) and after a CTL thinning operation was accomplished (rut_T) on newly assigned machine operating trails (Figure 1). An unequal variance t-test was used to compare between 4 values of rut_H within predicted sensitive areas, where the depth-to-water (DTW) index was below 1 m and 43 values, measured within areas with predicted high trafficability (see Section 2.4.1).

		rutH (cm)		rutT (cm)		Test Statistics for rutH
DTW	n	Mean	SD	Mean	SD	t-Value	p-Value
<1 m	4	2.65	2.50	5.17	1.54
>1 m	43	3.71	1.56	6.43	2.14	1.35	0.18

and Figure 1), usually assumed to indicate areas with a low bearing capacity and high susceptibility to soil displacement. Within this study, rut_H measured on these transects averaged at 2.65 cm. Besides, beneficial operational conditions were predicted for 43 measuring transects, where a similar extent of rut depth was estimated (Table 4). The relatively dry conditions, from the perspective of forest operations, might have caused a missing relationship between rut_H and DTW values (Figure 9A). In addition, it might impede a correlation between values of SMC_VOL and DTW values (p = 0.336, r_P = −0.14).

Similar to the observations made regarding DTW, the TWI was also not able to predict occurring rutting, since no correlation between log-transformed TWI and rut_H could be approved (Figure 9B).

4. Discussion

Within this study, a regularly scheduled forest operation was analyzed in two broadleaved dominated mixed stands, with machine traffic on 19 newly established permanent machine operating trails. Prior to the operation, measuring transects have been established along the machine operating trails. On each of these, terramechanical parameters were quantified in advance of the operation, using different equipment. Additionally, two cartographic indices were calculated. In order to validate the compiled parameters and indices, correlations with an occurring rut depth, measured after the first and final machine, were performed.

4.1. Soil Impact

Numerous factors, like wheel slippage, high traffic and the relationship between soil bearing capacity and a machine’s ground pressure, can lead to severe soil impacts, harming physical soil properties [58]. The long-lasting [59] impacts, disturbing the integrity of the soil, can be caused by all vehicles. Due to a compression stress and plastic deformation [60], soil impacts and displacement emerge.

After the off-road traffic by heavy forest machines (Table 1), values of soil compaction within this study averaged at 10%, related to the initial soil bulk density of 1.2 g cm⁻³. This is in line with the results reported in the meta-analysis by Ampoorter et al. [13], revealing an average increase in soil bulk density of 15% under similar operational conditions. Kozlowski [12] reported a critical soil bulk density of 1.4–1.6 g cm⁻³ to hamper root penetration into fine-textured soil [61]. It can be assumed that, in particular, the first machine passes prevailingly lead to soil compaction [13,62,63]. Consequently, additional machine passes result in a lower increase in soil bulk density, since preceding passes pre-compacted the soil already. In agreement, we have clear evidence for a higher soil compaction on measuring transects with a low initial soil bulk density (Figure 6B). Although soil compaction by forest machines might decrease with additional traffic frequency, rut depth increments can be enhanced after a certain compaction has been reached [64,65]. Soils are susceptible to compression, especially under moist and wet soil conditions [15,66], since the particle-to-particle bonding is low then [67]. The measured initial soil moisture content in the study area of only 29 vol% can be associated with a smaller extent of soil displacement [23,68] and might have impeded a vast soil displacement [64]. In this respect, Poltorak et al. [33] resumed deep ruts and a high soil compaction was more likely to occur in soils with a moisture content of more than 50 vol%, much higher than our measured 29 vol%.

Subsequently, ruts measured within this study appeared at a moderate extent, averaging at 6.3 ± 2.1 cm, with a maximum depth of 11.6 cm. Overall, only on 2 out of the 47 measuring transects, ruts deeper than 10 cm were observed. These results are in line with findings from other studies, reporting moderate rut depths between 3 cm [69], over 9–12 cm [8,63,70], and up to 20 cm [9,15,71], when measured on designated machine operating trails. In general, rutting that appeared within this study was in agreement with the recommendations of Owende et al. [14], who proposed to limit rut depths to a maximum of 10 cm to ensure an eco-efficient wood harvesting in Europe. Considering that the presented field study was conducted on newly assigned machine operating trails, the experienced rut formation was further below the internal guidelines of the state forest enterprise of Lower Saxony, stating a maximum tolerance value of a 15 cm rut depth on 90% of length of the machine operating trail [72]. These guidelines are set to ensure the technical trafficability of permanent machine operating trails for future operations, rather than actively protecting forest soils. Yet, since maintaining that the technical trafficability is inevitable to avoid the conversion of further stand areas to new machine operating trails, counteracting the soil conservation thoughts behind the permanent machine operating trail concept, it is fundamental to keep rut formation there as low as possible too [13,60]. However, it should be noted that the moderate extent of ruts experienced in this study after the crossing of a 24.5 Mg harvester, and subsequent timber extraction by a forwarder, is most likely supported by the relatively low logging volume of just 36 m³ ha⁻¹. This, in return, led to a low traffic frequency (Figure 1 and Figure 7B) with limited loads. It should be mentioned that rutting at the measuring transects was favored due to the conductance of this study, since operators were advised not to create brush mats in the transect areas, contrary to common recommendations to reduce rut formation [73].

4.2. Validation of Surveyed Tools

In general, observed rutting and soil compaction are in line with earlier findings [13,23]. Still, rut formations arising after forwarding possess an unfortunate nuisance within the current study, which is also reflected by the relatively weak correlation between rut_H and rut_F. The latter mentioned was probably deranged: (i) by changing moisture conditions, as rain fall occurred in the time gap between harvester and forwarder trafficked the study site; (ii) through the study’s arrangement within a real harvesting operation, where numerous measuring transects became unusable, since the wooden pegs, used to position the reference beam, were displaced, reducing our number of observations to 47 (from initial 90); and (iii) because loads and consequently total mass of the machine could not be captured during forwarding, reducing the potential independent variables in the linear models. Hence, linear models of the parameters to be validated were fitted with rut_H, created by the harvester, in order to account for the uncertainties (i–iii) within the validation of the surveyed approaches. Deviations from a standardized study design were assumed to be lower when using rut_H, compared to rut_T. Although, we assumed that one machine pass by the harvester (driving without GNSS-tracking) and short back and forth movements of the harvester may have led to multiple crossings on a measuring transect during the felling and processing, which is difficult to quantify at a given transect.

Although Sirén et al. [18,28] reported a Cone Index (CI), measured by means of a handheld Penetrologger, to be correlated with rut depth after a harvester pass (r_P = −0.27 [18] or r_P = −0.24 [28]), we were not able to statistically confirm this relationship within the current study (r_P = −0.12, p = 0.415). This instrument possesses a pivotal nuisance, as the measured values are highly dependent on user handling [21], which might have perturbed the correlation within this study. In addition, patterns of the upmost few centimeters of soil might not have been associated with actual bearing capacity, as organic matter content [74] or divergent soil moisture conditions can highly affect measured penetration resistance there. To avoid these uncertainties, we calculated a modified Cone Index, CI_MOD, which includes penetration resistance of between 10 cm and 20 cm soil depth, as compared to CI, considering the values measured in a depth between 0 and 15 cm. As shown in Figure 10, CI_MOD inversely responded to rut_H. Contrarily, common CI turned out to be the least accurate terramechanical parameter to predict rut_H, as summarized in Figure 8. Interestingly, the much easier to measure penetration depth (PD), herein also derived from the Penetrologger, turned out to be in a clear response with rut_H. A lower PD might be driven by increased coarse fraction [75] or a more intense root network within the soil, which both can be considered to impede rutting.

Besides penetration depth, the dual-mass dynamic cone penetrometer also uses a metric scale (DCP). Measurements of this relatively simple mechanical device appeared to show robust predictions of rut_H and resulted in a more reliable r_P (0.11), with total rut depth, rut_T, compared to CI. The latter-mentioned, curiously, showed a positive trend to rut_T, meaning that deeper ruts tended to occur where highest values of CI were captured (p = 0.15). During the application of the dual-mass dynamic cone penetrometer, the hammer is raised carefully without lifting the cone before it is allowed to drop [76]. By doing so, an artificial bias during handling can be avoided almost certainly. We assume that this, in combination with a larger cone, resulted in a better prediction of rut_H, compared to the Penetrologger. However, it should be noted that the penetrometer is an instrument with a weight of 12.8 kg, limiting its application for individual users at extended stand levels.

Although seldomly used in research, shear strength (τ), measured by the use of a vane meter, led to rather reliable results, regarding the correlation of τ and rut_H. Thereby, r_P of −0.25 exceeded the coefficients achieved by CI. This correlation is in close agreement with the findings of Jacke et al. [77], who promoted the usage of a vane meter to measure shear strength and to assess temporal soil trafficability, as already verified [30].

In summary, the surveyed terramechanical parameters PD, DCP and τ (in decreasing order according to revealed r_P) bear a high potential towards the prediction of rutting.

4.3. Correlations with Cartographic Indices and Soil Moisture

Nonetheless, these measurements are time-consuming and might not be doable in day-to-day work by harvesting operators and entrepreneurs. Subsequently, remotely analyzed trafficability maps, covering the whole area of operation, could match the demands of practitioners, presuming that those maps show a valid prediction of expectable rutting. Both used cartographic indices simulate water saturation or soil wetness, as it would accumulate in depressions, sinks or if a sufficiently sized area drains into a spot. However, measured spatial distribution of soil moisture content (SMC_VOL) could not be reflected by depth-to-water (DTW) maps within the present study, as shown in Figure 5. This might have also caused the contradictory response between DTW index and observed rut_H, which occurred completely independent of DTW (Figure 9A). The analyzed DTW maps, which were calculated with a flow initiation area of 4 ha [36,37], would represent a low overall moisture [29]. The created index is usually converted to a dummy variable, with a threshold of 1 m [51], where values below this threshold are assumed to predict wet and sensitive areas. By doing so, four measuring transects would have been assessed to have an inferior bearing capacity due to the assumed water saturation. Measured rut_H and rut_T at these four transects showed low means of 2.7 cm or 5.2 cm, respectively, similar to values measured on the remaining transects (Table 4). However, considering the generally acceptable extents of rut_T [14,72], the DTW index reasonably assessed actual operational conditions for 43 out of 47 transects, as favorable trafficability was predicted there. In addition, the occurring precipitation between felling and forwarding may limit the meaningfulness of the made observations. The validity of this approach could be restricted, since soil moisture increased particularly at dryer transects (Figure 6A), thereby smoothing the spatial distribution of moisture. Hence, the prediction of wet or water-saturated areas might be jeopardized. Along with the DTW, also the second surveyed cartographic index, topographic wetness index, revealed no reasonable correlation with rut_H (Figure 9B).

The inferior performance of both surveyed indices might be driven by two main aspects: (i) Due to the systematic allocation of measuring transects, some of these were positioned at former machine-operating trails, which could be observed at both forest stands, although they were assumed to be un-trafficked. Whereas terramechanical test procedures capture the divergent soil conditions there, trafficability maps are not able to include these effects. (ii) The influence of soil moisture on rutting has been reported extensively [23,33], with inevitable influence on soil displacement. However, within this study, the measured rutting turned out to be independent of soil moisture content (SMC_VOL), which is simulated by DTW and TWI in the first place. The relatively low SMC_VOL (from the perspective of forest operations) during the field trials, with only six observations above 40 vol%, could have limited the effect of SMC_VOL [74] on the appearing rut formation and consequently the predictability of rut depth by means of both cartographic indices.

4.4. Prediction of Rutting

As none of the surveyed tools were able to predict final rutting to a satisfying level at the given conditions, alternative approaches are requested to match the desired prediction of soil impact [20,21]. Findings of Sirén et al. [18,28] approved the rut formation of the first machine pass, and in our case, by the harvester (rut_H), to predict further rutting caused by the consecutive passes of the forwarder. In agreement with these findings, rut_H of the present study turned out to be the best indicator for total rut depth after 2 to 5 machine passes, possessing a highly significant correlation among both variables (Figure 7A). Under the given conditions, higher traffic frequency did not lead to additional rut depth increments (Figure 7B), which supports the overall assumptions of this survey. The results give a further impetus to use rut_H with regard to the prediction of rutting. In the modelling of Sirén et al. [18], variables like SMC and CI did not further increase accuracy when rut_H had been included already. In agreement with this, rut_H acted as a pivotal influence, explaining the highest share of deviations in the present study, where additional variables did not increase the model’s quality (not shown). An exchange of information between both operating machines could contribute to the optimal routing of the subsequently driving forwarder and thereby support an eco-efficient forestry through the avoidance of areas with a high susceptibility to rut formation. Making it useful for operational purposes, Lidar-based sensors could detect actual rutting, created by an operating harvester, evolving the real-time generation of trafficability maps as a part of normal forestry operations [26]. Another sensor-based solution aims at rolling resistance, which has been recently assumed to be able to assess machine-operating trail trafficability [27]. In agreement, Salmivaara et al. [78] reported the high potential of CAN-bus-derived rolling resistance to be used for predictions of trafficability across complex landscapes and changing conditions.

These more operation-based monitoring approaches in particular have potential to support the management of extraction activities, consecutive to the harvester’s felling and processing, in order to prevent further soil impact whilst forwarding. Such technologies could support harvesting operations without fixed machine operating trails, as conducted in many plantations of the tropics or boreal forests. If the harvester detected areas with low bearing capacity, timber extraction could be conducted on alternative routes, in order to avoid severe soil impact. However, this does not solve the problem of the impact caused by the first machine pass, which requires a reliable pre-operation assessment of the trafficability situation at both, permanent machine-operating trails with pre-compaction, and previously non-trafficked areas.

5. Conclusions

The findings of the current study can be used to answer the posed research questions. One of these questions addressed the potential prediction of total rut depth after CTL operations, based on rutting after an initial machine pass. The results give clear evidence that rut depth caused by the first machine pass of the harvester remains a reliable indicator of further rut formation through up to four consecutive machine passes in current site conditions, and should therefore be cautiously monitored during the felling. However, for operational planning purposes, assessment of the logging sites prior to machine traffic is mandatory to ensure efficient machine utilization and to avoid negative impacts on soils or the technical trafficability of permanent machine operating trails in the first way. In this respect, the auspicious remote sensing-based cartographic indices DTW and TWI seem promising towards an operational usage, presuming a high accuracy of rut prediction. However, both used indices did not prove to be reliable alternatives compared to conventional terramechanical in-situ measurements of actual site conditions, at least during the experienced situation with already long-lasting precipitation deficits. The surveyed terramechanical test procedures possess specific advantages and drawbacks, but have not been compared with each other so far. This begs the question of the tool-specific ability to assess operational conditions in terms of trafficability, and in turn the ability to predict ruts. Surprisingly, the frequently used site assessment indicator CI showed a low certainty of rut depth prediction, compared to the less commonly applied dual-mass dynamic cone penetrometer and vane tester. Finally, none of the investigated methods can be fully approved or disapproved due to the limited scale of observations and the intensity of monitored machine traffic in this study. However, the results clearly revealed that the various approaches, already available for operational planning, would lead to a different assessment of the situation. Thus, state-of-the-art applications like cartographic indices need further adaptation to increase the predictive power on the one hand, but also parameters measured by more simple instruments such as the DCP, should be further considered as an alternative, or at least as a supplementary aid to determine current site conditions and risks of trafficability.