Impact of Skidding and Slope on Grapple Skidder Productivity and Costs: A Monte Carlo Simulation in Eucalyptus Plantations

Abstract

Background: In the context of mechanized timber harvesting, alterations in technical parameters, such as skidding distance and terrain slope, have the potential to influence the productivity and production costs associated with the self-propelled grapple skidder. Furthermore, these variables are inherently uncertain, which could potentially cause forest managers to make inaccurate decisions. The objective was to analyze whether four skidding distances and two slope classes influence the productivity and production costs of the grapple skidder in Eucalyptus planted forests from a stochastic perspective using the Monte Carlo method. Methods: Productivity was estimated using the time study protocol. To calculate the cost per scheduled hour of the grapple skidder, both fixed and variable costs were considered, and subsequently, the production cost was determined. Results: The mean productivity of the grapple skidder on flat slopes was 114.35 m³ h⁻¹, while on wavy to strong wavy slopes it was 80.43 m³ h⁻¹. In flat slopes, considering all skid distance ranges, the mean production cost was 0.82 USD m⁻³, while in wavy to strong wavy slopes it was 1.48 USD m⁻³. The mean values for operator labor costs and fuel account for 58.1% of the cost per scheduled hour of the grapple skidder. Conclusions: The mean productivity of the grapple skidder in Eucalyptus planted forests decreased with increasing skidding distance in both slope classes but was 29.7% lower on wavy to strong wavy slopes compared to flat slopes. The mean production cost of the grapple skidder during timber skidding on flat slopes is 80.0% lower than on wavy to strong wavy slopes. For future investigations, the impact of other slope classes, skid distances, and silvicultural aspects on productivity and production costs can be considered from a stochastic perspective using the Monte Carlo method.

1. Introduction

In the context of mechanized timber harvesting in planted forests, the efficiency and production costs associated with the use of self-propelled grapple skidders for skidding timber bundles can vary depending on various operational factors, such as the slope of the terrain and the distance covered during the skidding process. Thus, utilizing stochastic methods can enable forest managers to comprehend the impact of variations in these factors on the key performance indicator. These methods are known for their ability to provide reliable data, facilitate autonomy in strategic planning, and support informed decision-making.

Planted forests currently cover 294 million hectares worldwide, representing 7% of total forest area. The area has grown annually by just under 1% from 2015 to 2020 [1]. In Brazil, forest plantations cover an extensive total area of 9.5 million hectares, whereby the majority, 76.8%, accounts for Eucalyptus spp., resulting in a total of 7.3 million hectares. Eucalyptus planted forests are usually harvested using mechanized systems, which function as raw material sources for Brazil’s forest-based industries [2].

The management of mechanized timber harvesting activities entails selecting a system for adoption by forest-based companies that relies heavily on factors such as soil and climate conditions, forest traits, assortments, skilled labor, and financial availability. The full-tree harvesting system involves felling trees without their roots and creating bundles of timber that are then extracted from the plot and transported to the forest roadside for subsequent sectioning [3,4,5].

When mechanized methods are used in the full tree harvesting system, the grapple skidder is a frequently employed machine in the logging industry, serving as a self-propelled, articulated forestry vehicle with a grapple-style attachment for skidding timber bundles and employing pneumatic wheels. In forestry management and planning, the primary factors for appraising mechanized timber harvesting operations are machine productivity and production costs [6,7,8].

The productivity of a grapple skidder may vary depending on several factors, such as the machine’s model, operator’s level of experience, terrain slope, average tree diameter and volume, environmental conditions, and skidding distance. Production costs, however, may be influenced by the residual and purchase value of the grapple skidder, its useful life, fees, interest and taxes, fuel consumption, maintenance, and mechanical availability [9,10,11]. Predicting the productivity and production costs of grapple skidders is uncertain due to their constant variations, which could potentially cause forest managers to make inaccurate decisions.

The act of disregarding uncertainty when making decisions is predicated on the assumption that the future value of variables will align with their expected or mean value. However, it is uncommon for variables to align perfectly with their expected value, leading to suboptimal harvest decisions during the early planning horizon and hindering the achievement of planning objectives. Considering the scope of deterministic models, it is important to note that the solution bounds of the stochastic model are not predictive of the future. However, they do provide a range of possibilities for further fits [12,13,14].

The Monte Carlo method, also referred to as the statistical simulation approach, is a widely suggested technique for analyzing uncertainty. This method introduces some level of randomness into the model being used, with each simulation (i.e., scenario of uncertainty) generating a distinct evaluation value for each criterion using the corresponding probability distribution. This approach represents a novelty in the scientific literature, in contrast to deterministic methods often adopted in timber harvesting studies, because it allows the evaluation of the stability of the results by introducing a stochastic element in the input variables. The introduction of this element emulates the innate uncertainty of the model’s inputs, thereby supporting decision-makers in making more precise decisions [15,16,17].

Despite the increasing significance of conducting uncertainty analysis in forest operations based on scientific research, numerous forestry sector managers continue to disregard it. This study aims to assess the impact of slope classes and skidding distances on both the productivity and production costs of grapple skidders within Eucalyptus planted forests through a stochastic perspective utilizing the Monte Carlo method.

The article is organized as follows: First, the theoretical background is presented to show that the prediction of productivity and production costs of grapple skidders could be improved if presented in a stochastic framework using the Monte Carlo method. Then, the survey data and the methods used to estimate the stochastic model are presented. Successively, the estimated productivity and production costs and their probability distributions are presented. Finally, the main findings, limitations, and possible extensions of the approach are discussed.

2. Materials and Methods

2.1. Study Area

The grapple skidding activity was conducted in a 72-month-old Eucalyptus platyphylla planted forest with an average individual volume of 0.24 m³ and 3.0 m × 2.0 m spacing. The forest is located at geographical coordinates 23°06′ S and 48°36′ W in the State of São Paulo, Brazil. In accordance with the Köppen-Geiger classification system, ref. [18] has identified the region’s climate as humid subtropical (Cwa), which is typified by arid winters and hot summers. The region experiences an average annual rainfall of 1350 mm, an average relative humidity of 64.0%, and an average temperature of 20.0 °C. The soil was classified as dystrophic Red Yellow Latosol, medium texture [19].

2.2. Experimental Details

The mechanized harvesting of timber was conducted using the whole-tree system under conditions of shallow management. The trees were felled using a self-propelled forestry machine, defined as a feller-buncher by the International Organization for Standardization [20], which placed two bunchers, containing an average of 20 trees each, on the ground for subsequent skidding.

The skidding of timber entailed the grouping of four bunchers from the felling site, which were then placed at an angle of 90° to the forest roads by a self-propelled forestry machine, designated a grapple skidder by the International Organization for Standardization. Subsequently, the bunchers were processed by the self-propelled grapple saw with the operating standard outlined in the International Organization for Standardization [21].

The grapple skidder (John Deere 848H, Deere and Company, Moline, IL, USA) had a nominal power of 149 kW and was equipped with pneumatic wheels measuring 35.5 mm in width. The grapple has a radius of 32 inches and operates at 45 bar pressure. It is equipped with a hydraulic grapple with a capacity of 1.5 m², weighs approximately 17,826 kg, and has accumulated 15,657 h of use. For the purposes of this study, only one operator was considered. This individual was male, 27 years of age, and had 24 months’ experience in skidding timber with the aforementioned forestry machine.

Following guidelines from the scientific literature [22,23,24], timber was removed from the felling site in four skidding distances (SD), which were classified as follows: 0–50 m; 51–100 m; 101–150 m; 151–200 m. The two slope classes (SC) of the study area were determined in the field using clinometers and were classified by [19] as SC 1 (flat, ≤3.0%) and SC 2 (wavy to strong wavy, 8.0–45.0%).

2.3. Time Study Application

The sizing of production resources utilized in forestry operations is traditionally conducted through the implementation of the time study protocol. As outlined by [25], the protocol was implemented through the continuous method, employing a manual digital stopwatch with one-second accuracy and without interruption.

Previously, reference [26] conducted a pilot study for 372 min to observe the machine elements (ME) that constituted the operational cycle. These were identified as empty travel (ET), buncher loading (BL), travel loaded (TL), and buncher unloading (BU).

The sample size (Equation (1)) for the slope classes SC 1 and SC 2 with a confidence level of 95.0% was determined from the timing of a random sample taken preliminarily, based on, $z_{\alpha /2}^2$ which is the critical value of the normal distribution at $\alpha /2$ (e.g., for a confidence level of 95%, α is 0.05 and the critical value is 1.96), σ², which is the sample variance obtained in the pilot study, and E, which is the margin of error, that is, the desired level of precision [27,28]. Sample sufficiency calculations have been used in recent studies in the field of mechanized timber harvesting [29,30].

$$n={\displaystyle \frac{z_{\alpha /2}^2{\sigma }^2}{E^2}}$$

(1)

where $n$ is the sample size, $z_{\alpha /2}^2$ is the tabulated value of z for standard normal distribution, ${\sigma }^2$ is the sample variance, and $E^2$ is the margin of error of the estimate.

2.4. Productivity Study

The productivity of the grapple skidder (Equation (2)) was obtained in accordance with the methodology proposed by [31], whereby the volume of skidded timber and the effective time of the operational cycle were recorded.

$$P={\displaystyle \frac{v}{t}}$$

(2)

where $P$ is the productivity per effective hour (m³ h⁻¹), $v$ is the volume of skidded timber (m³), and $t$ is the effective time (h).

2.5. Economic Management

The cost per scheduled hour of the grapple skidder (USD h⁻¹) was calculated based on the machine’s scheduled hour, inclusive of all times, and the utilization rate (

Table 1. Cost assumptions and hourly rates for the John Deere 848H grapple skidder.

Factor	Unit	Value
Initial investment	USD	296,330.68
Residual value	USD	59,266.14
Fuel consumption	L h−1	25.92
Fuel price	USD L	1.24
Rated motor power	kW	149
Economic life	h	30,000
Number of days worked per year	d	283
Number of shifts per day	d	3
Scheduled hours per shift	h	8
Utilization rate	%	74.0
Estimated service life of the tire set	h	5000
Operator’s basic salary	USD h−1	13.93
Social charges	%	134.0

), in accordance with the methodology proposed by [32]. In this costing method, fixed costs were defined as those that did not depend on the use of the grapple skidder. These included depreciation, return on capital applied to the acquisition of fixed assets, insurance, shelter, property tax, and transportation of the grapple skidder. As variable costs, those components associated with the intensity of grapple skidder use were weighted. This included the costs of fuel, lubricating oils and greases, maintenance and repairs, tires, operator labor, and 5.0% overhead, which was calculated from the fixed costs.

Given the financial contribution from various sources of capital for the acquisition of the grapple skidder, the opportunity cost rate for remunerating this capital was calculated using the Weighted Average Cost of Capital (WACC). This method reflects the participation of each source in a weighted average of the marginal cost of capital, as described by [33] and illustrated in Equation (3):

$$WACC={\displaystyle \frac{\left[D k_D \left(1-CT\right)+E k_E\right]}{\left(D+E\right)}}$$

(3)

where $D$ is the creditor’s cost of capital, $k_D$ is the interest rate on the creditor’s capital, $CT$ is the corporate tax rate, $E$ is the market value of the debt, and $k_E$ is the interest rate on equity capital.

As the interest rate on equity is not directly observable in the financial market, it was calculated by inference from the Capital Asset Pricing Model (CAPM). As stated by [34], the country risk premium was incorporated into the analysis, given that the study was conducted in an emerging economy (Equation (4)):

$$k_E=R_f+\beta l\left(R_m-R_f\right)-{\mathsf{\Omega }}_{Br}$$

(4)

where $R_f$ is the rate of return on a risk-free asset, $\beta l$ is the asset’s systematic risk coefficient, $R_m$ is the expected rate of return for the market portfolio, $R_m-R_f$ is the market risk premium, and ${\mathsf{\Omega }}_{Br}$ is the country risk premium.

In order to ascertain the systematic risk coefficient of the asset in question, the average unlevered beta of publicly traded forestry companies in Brazil was employed as a point of reference. The following companies were considered: Companhia Melhoramentos (São Paulo, Brazil), Dexco S.A. (São Paulo, Brazil), Eucatex S.A. Indústria e Comércio (São Paulo, Brazil), Klabin S.A. (São Paulo, Brazil), and Suzano Papel e Celulose S.A. (São Paulo, Brazil) [35]. The proportion of assets financed by debt and the corporate tax rate of 34.0% were considered, as recommended by [36], to allow for the capture of tax benefits resulting from interest payments. The real levered beta was determined to be 0.77, with the proportion of assets financed by debt assumed to be 44.5%.

A premium of 3.3% has been appended to the cost of capital for creditors in countries with a speculative credit rating of Ba2, according to [37]. The risk-free interest rate of 5.1% was defined using the geometric mean of the period between 1 February 1962 and 29 December 2023 of the annual return on treasury bonds with a 10-year maturity rate provided by the [38]. The country risk premium of 3.8% was calculated using the geometric mean of the historical series of Brazil risk between 29 April 1994 and 1 December 2023, in conjunction with the Emerging Markets Bond Index Plus, as published by [39].

The anticipated rate of return for the market portfolio was 4.9%, as indicated by the S and P Global Timber and Forestry Index published by [40]. A market risk premium of 0.2% was calculated by employing the risk-free interest rate and the anticipated rate of return for the market portfolio. The cost of equity was determined to be 8.4%. By incorporating the proportion of third-party capital and the proportion of own capital, which was 55.5%, it was feasible to ascertain the opportunity cost rate of 7.3%.

The production cost of the grapple skidder was calculated using the ratio between the cost per scheduled hour and productivity [30], or, in other words, the volume of timber skidding per effective hour of work (Equation (5)):

$$C_{pm}={\displaystyle \frac{C_{hp}}{P}}$$

(5)

where $C_{pm}$ is the production cost of the grapple skidder (USD m⁻³) and $C_{hp}$ is the cost per scheduled hour of the grapple skidder (USD h⁻¹).

2.6. Stochastic Modeling

The stochastic model was developed based on stochastic (random) variables, or in other words, the inputs of the mathematical models, namely the times of the machine elements and the volumes of the skidded timber, to which probability distributions were fitted to the sample sets in order to determine the probability distributions that best described the data, i.e., to ensure the adequacy and validation of these data using the Bayesian Information Criterion (BIC), which is widely used to evaluate the quality of fit of a given model. Therefore, according to [41,42], the BIC was calculated based on the logarithm of the likelihood function. Nevertheless, we described the uncertainties using probability distributions. Furthermore, the triangular probability distribution was assigned to the base value of the fixed cost components with a variation of ±15.0%, according to the recommendations of [43].

A Monte Carlo simulation was conducted to obtain a range of values using the @Risk version 8.8.1 software [44], with 100,000 pseudo-random numbers generated [45]. The Mersenne twister pseudorandom number generator was employed in the simulation process, in accordance with the methodology proposed by [46]. The initial parameters for the mathematical model were set [47]. The Kolmogorov-Smirnov (K-S) test was applied at a significance level of 1.0% to assess the normal probability distribution of the data [48,49].

The outputs of the mathematical models were the cost per scheduled hour, productivity, and production cost of the grapple skidder [50]. As the data were non-parametric, Spearman’s rank correlation coefficient was calculated, with the parameter represented by $p_s$, at a significance level of 5.0% [51,52]. In order to interpret the monotonic relationships between the machine elements and the volume of skidded timber, the intensity of association between inputs and outputs was based on the methodology proposed by [53].

3. Results and Discussion

3.1. Stochastic Analysis of Machine Elements

In determining the minimum sample size, 86 and 275 cycles were considered for slope classes SC 1 and SC 2, respectively. Consequently, the number of operational cycles observed exceeded the minimum sample size of 76 cycles for both slope classes. The total effective working time was 22 h and 6 min, allowing 1673.5 m³ of timber to be skidded.

The time study protocol has become an important tool for forest resource assessment because it allows quantification of the time required to perform each activity. In addition, when combined with sample sufficiency, it provided statistical reliability. The observation of a greater number of cycles than the minimum required made it possible to reduce the margin of error of the estimate to 4.0%, thus increasing the credibility of the results, as discussed by [54].

The empirical conditions for skidding timber in SC 1 yielded an average operating cycle time of 4 min and 45 s, while in SC 2, this time was 4 min and 42 s. The TL machine element had the highest average representation of the total time (39.3%), followed by the ET (34.5%), BL (17.8%), and BU (8.4%) elements. These results are consistent with the conclusions of [55], which indicated that displacement elements required more time than the loading and unloading of bundles.

In the SC 1 slope class, the best fits of the normal distribution were observed for the ET machine element in the 0–50 m skid distance range, TL in the 51–100 m range, and BL in the 151–200 m range. As noted by [56,57], this probability distribution is significant because it represents real-value random variables whose distributions are not known.

The uniform probability distribution showed the best fit for the machine element BU in the 0–50 m skid distance range and ET in the 51–100 m and 101–150 m ranges. As postulated by [58], this type of fit is more discriminating with parameters for the lower and upper limits. Consequently, the times of these machine elements were based on these specified limits, i.e., probable values.

Considering the triangular probability distribution, the optimal fits for the machine element TL were observed in the ranges of 0–50 m, 101–150 m, and 151–200 m. According to [59,60], this distribution is used because of its simple geometric shape when there are uncertain parameters in the process. Consequently, it is plausible that this shape is due to the heterogeneity of the forestry activity, which includes distinctive characteristics from the formation of bunchers to the skidding at the ends of stands.

The exponential probability distribution has a wide range of applications in various fields [61,62]. It was observed that the proposed fit showed a better fit to the BL machine element data in the 0–50 m skid distance range and to the BU in the 51–100 m, 101–150 m, and 151–200 m ranges. Consequently, the application of this adjustment resulted in improved accuracy in the mathematical modeling of this machine element.

The machine element BL, located within the skid distance range of 51–100 m, exhibited a probability distribution that was consistent with the Laplace distribution. As stated in [63,64], this is one of the earliest distributions in probability theory and is used when the scale parameter is known and its finite sample exhibits asymptotic properties. In light of the above considerations, given the finite number of observations, it was possible to assign a probability distribution that closely approximates a normal distribution to the time required to load the bunchers.

When this machine element was analyzed within the 101–150 m skid distance range, it was found that the logistic distribution provided the most accurate fit, thus supporting the attribution. This finding is consistent with the conclusions of [65,66], which demonstrated the effectiveness of this distribution in the presence of random parameter uncertainties. Accordingly, the implementation of the time study protocol may inherently involve certain uncertainties regarding the time series of observations and the operational factors involved.

Considering the four skid distance ranges, the machine elements followed probability distributions with specific behaviors (

Table 2. Probability distribution and descriptive statistics of the times in minutes of the grapple skidder machine elements in the flat slope class as a function of the skidding distance ranges.

SD 0–50 m	ME	Distribuition	Minimum	Maximum	Mean	Modal	S.D.	Percentile 5	Percentile 95	BIC	K-S
	ET	Normal	−∞	+∞	0.6204	0.6204	0.2237	0.2524	0.9884	1.9576	0.0835
	TL	Triangular	0.2333	1.8690	0.7786	0.2333	0.3855	0.2747	1.5032	24.0973	0.1264
	BL	Exponential	0.0795	+∞	0.5097	0.0795	0.4302	0.1016	1.3682	14.6578	0.1158
	BU	Uniform	0.0614	0.5219	0.2917	-	0.1329	0.0844	0.4989	−25.0239	0.3300
SD 51–100 m	ME	Distribuition	Minimum	Maximum	Mean	Modal	S.D.	Percentile 5	Percentile 95	BIC	K-S
	ET	Uniform	0.5306	1.9028	1.2167	-	0.3961	0.5992	1.8342	22.2593	0.1371
	TL	Normal	−∞	+∞	1.2278	1.2278	0.2838	0.7610	1.6945	13.0042	0.0729
	BL	Laplace	−∞	+∞	0.6000	0.6000	0.2511	0.1912	1.0088	4.6411	0.1164
	BU	Exponential	0.0753	+∞	0.2594	0.0753	0.1841	0.0848	0.6267	−23.5842	0.3206
SD 101–150 m	ME	Distribuition	Minimum	Maximum	Mean	Modal	S.D.	Percentile 5	Percentile 95	BIC	K-S
	ET	Uniform	1.0897	2.7769	1.9333	-	0.4870	1.1741	2.6926	34.8368	0.2185
	TL	Triangular	1.1333	3.0825	1.7831	1.1333	0.4594	1.1827	2.6467	37.6380	0.1478
	BL	Logistic	−∞	+∞	0.6927	0.6927	0.1803	0.4000	0.9853	−7.4970	0.0878
	BU	Exponential	0.0785	+∞	0.2094	0.0785	0.1309	0.0852	0.4705	−47.2225	0.2655
SD 151–200 m	ME	Distribuition	Minimum	Maximum	Mean	Modal	S.D.	Percentile 5	Percentile 95	BIC	K-S
	ET	Triangular	1.7667	2.7876	2.1070	1.7667	0.2406	1.7925	2.5593	3.8946	0.2176
	TL	Triangular	1.4670	2.7000	2.2890	2.7000	0.2906	1.7427	2.6688	7.4209	0.1465
	BL	Normal	−∞	+∞	0.7852	0.7852	0.3021	0.2882	1.2821	7.3910	0.1829
	BU	Exponential	0.0891	+∞	0.1653	0.0891	0.0762	0.0930	0.3174	−16.1514	0.4383

ME = machine elements, ET = empty travel, TL = travel loaded, BL = buncher loading, and BU = buncher unloading.

), namely Exponential corresponded to 25.0% of the probability distributions, triangular (25.0%), normal (18.8%), uniform (18.8%), Laplace (6.2%), and logistic (6.2%). The exponential distribution allowed the identification of the machine elements with values close to zero, which therefore required the least time of all the others.

The triangular probability distribution allows forecasting without the necessity of assuming symmetry around the mean value of the element times. In addition, reference [67] states that both the exponential and triangular distributions are more commonly used because they correspond to the cognitive processes and decision-making styles of the majority of managers and are easy to implement.

In the context of the SC 2 slope class, the ET machine element exhibited a superior distribution fit in the 0–50 m skid distance range, as evidenced by the Weibull distribution. This distribution, as postulated by [68,69], is often used to model phenomena that exhibit monotonous failure rates. In other words, the time required to move the grapple skidder without dragging tree bundles over distances of less than 50 m was found to be constant.

However, within the 151–200 m skid distance range, the logistic distribution proved to be the optimal fit for this machine element. As stated by [70], this probability distribution has become a standard tool for investigating the probability that an event can be influenced by one or more explanatory variables. Consequently, the only explanatory variable was the empty travel speed, which was found to be higher than the loaded travel speed.

The gamma distribution was found to be an optimal fit for the TL machine element in the 0–50 m skid distance range. This range is widely recognized as a very common and adaptable domain, with applications in various scientific disciplines. In light of the above evidence, the model was shown to be reliable when applied to the data series of travel times with tree bundles along the shortest skidding corridor.

The exponential distribution was found to be the best fit for the BL and BU machine elements in the 0–50 m, 51–100 m, and 151–200 m skid distance ranges. As described in [71,72], this probability distribution is used to describe time. This allows for the modeling of continuous processes, which is the chosen time method.

Of the five distributions that showed optimal fit to the wavy to strong wavy slope data, the triangular probability distribution showed the most accurate fit to the ET, TL, and BU machine elements in the 51–100 m and 101–150 m skid distance ranges, and to the TL and BL machine elements in the 151–200 m skid distance range. This distribution was used by [73,74] in cases where the correlation between the variables is known but the data is uncertain. Consequently, one of the uncertainties was the value of the standard deviation of the measurement times of these machine elements, which is correlated with the estimates, mainly between the number of trees loaded and skidded.

In the four skid distances for SC 2, the machine elements resulted in different probability distributions (

Table 3. Probability distribution and descriptive statistics of the times in minutes of the grapple skidder machine elements in the wavy to strong wavy slope class as a function of the skidding distance ranges.

DA 0–50 m	ME	Distribuition	Minimum	Maximum	Mean	Modal	S.D.	Percentile 5	Percentile 95	BIC	K-S
	ET	Weibull	0.0992	+∞	0.6670	0.5424	0.3034	0.2390	1.2227	54.1309	0.0543
	TL	Gama	0.1154	+∞	0.5293	0.3792	0.2492	0.2199	1.0055	−7.9507	0.0586
	BL	Exponential	0.0642	+∞	0.3423	0.0642	0.2781	0.0784	0.8974	−51.2124	0.1130
	BU	Exponential	0.0477	+∞	0.3023	0.0477	0.2546	0.0608	0.8105	−70.9973	−0.1572
DA 51–100 m	ME	Distribuition	Minimum	Maximum	Mean	Modal	S.D.	Percentile 5	Percentile 95	BIC	K-S
	ET	Triangular	0.4857	1.5974	0.9777	0.8500	0.2314	0.6280	1.3935	1.5935	0.0755
	TL	Triangular	0.2910	3.1845	1.4363	0.8333	0.6279	0.5711	2.6013	13.8607	0.0712
	BL	Exponential	0.0935	+∞	0.5471	0.0935	0.4536	0.1168	1.4523	39.8126	0.1051
	BU	Triangular	0.0833	0.7126	0.2931	0.0833	0.1483	0.0993	0.5719	−67.8813	0.2165
DA 101–150 m	ME	Distribuition	Minimum	Maximum	Mean	Modal	S.D.	Percentile 5	Percentile 95	BIC	K-S
	ET	Triangular	0.6726	2.6391	1.4539	1.0500	0.4261	0.8652	2.2438	51.6372	0.0911
	TL	Triangular	0.6525	3.6176	2.0178	1.7833	0.6109	1.0620	3.0961	88.8248	0.1392
	BL	Exponential	0.0975	+∞	0.9388	0.0975	0.8413	0.1407	2.6178	82.3601	0.1655
	BU	Triangular	0.0833	0.9798	0.3821	0.0833	0.2113	0.1060	0.7793	−15.2818	0.1750
DA 151–200 m	ME	Distribuition	Minimum	Maximum	Mean	Modal	S.D.	Percentile 5	Percentile 95	BIC	K-S
	ET	Logistic	−∞	+∞	1.8272	1.8272	0.4864	1.0375	2.6168	49.1824	0.1259
	TL	Triangular	1.5333	4.3931	2.4866	1.5333	0.6741	1.6057	3.7536	58.8751	0.1248
	BL	Triangular	0.1500	2.4588	0.9196	0.1500	0.5442	0.2085	1.9425	46.8646	0.2463
	BU	Exponential	0.0765	+∞	0.2748	0.0765	0.1983	0.0867	0.6705	−27.1150	0.2675

ME = machine elements, ET = empty travel, TL = travel loaded, BL = buncher loading, and BU = buncher unloading.

): triangular corresponded to 50.0%, exponential (31.3%) and weibull, gamma, and logistic together accounted for 18.7%. Therefore, the triangular probability distribution, which only requires the minimum, most probable, and maximum values, can be used in the absence of data on the machine elements that make up the grapple skidder’s operating cycle, especially since it is a simple distribution, as described by [75].

3.2. Stochastic Analysis of Productivity

Productivity of self-propelled forest machines is a key performance indicator (KPI) [76]. Therefore, it is essential to consider the management of consumer units for both macro- and micro-level planning of timber harvesting activities.

It should be noted that the mean productivity of the grapple skidder decreased with increasing skidding distance in both slope classes. The mean productivity in slope class SC 1 was 114.35 m³ h⁻¹, while in SC 2 it was 80.43 m³ h⁻¹. Authors [77,78,79] have reported that among the factors influencing the productivity of the grapple skidder is the skidding distance. Therefore, the grapple skidder showed a reduction in productivity when it spent more time traversing longer distances, regardless of the slope class.

Another factor that has been identified as potentially influencing the productivity of grapple skidders is the slope of the terrain. This is a view that is supported by the findings of [80,81,82]. In the study by [83], the grapple skidder showed a 23.0% reduction in productivity in the steepest part of the plot (characterized by steep slopes) compared to the plot located on a wavy to strong wavy slope.

Consequently, the increased engine power required to traverse the terrain resulted in a reduction in the load capacity of the grapple skidder. This, in turn, led to a decrease in the mean productivity observed in slope class SC 2, with a mean reduction of 29.7%. This finding is corroborated by the studies of [84,85], which highlight the influence of terrain slope, engine power, drag load, and traffic conditions on the productivity of this machine.

The probability distribution that best fit the productivity data was the normal distribution (

Table 4. Mean productivity (m³ h⁻¹) and probability distribution of the grapple skidder as a function of slope class and skidding distance range.

Slope Classes	Skidding Distance Ranges (m)	Mean Productivity (m3 h−1)	Probability Distribution
SC 1	0–50	149.01 ± 63.62	Normal
	51–100	129.55 ± 29.34	Normal
	101–150	100.37 ± 22.43	Normal
	151–200	78.48 ± 16.37	Normal
SC 2	0–50	111.52 ± 55.69	Normal
	51–100	83.31 ± 19.28	Normal
	101–150	68.35 ± 27.18	Normal
	151–200	58.53 ± 21.92	Normal

). The result of this distribution was a symmetric curve, indicating that the grapple skidder productivity provided a model with a peak in the center and tails that tapered towards the mean of the data series obtained from the minimum sampling sufficiency.

As noted by [86,87], this distribution is considered one of the most important theoretical continuous probability distributions, applied to a wide variety of phenomena and used extensively for the theoretical advancement of statistical inference. In this context, the application of this distribution allowed the statistical validation of productivity.

When analyzing the monotonic relationship between skid volume of timber and productivity in the two slope classes and in all skid, distance ranges, the paired observations yielded positive correlations (p-value < 0.05), meaning that as skid volume of timber increased, productivity changed positively. Following [88,89], this correlation provided a result with greater robustness and reliability, as it is considered a non-parametric classification correlation, i.e., it does not depend on whether the data were linear or not.

Based on the above evidence, it can be concluded that in slope class SC 1, the correlations in question can be interpreted as moderate (${\mathrm{p}}_{\mathrm{s}}$ = 0.53) in the 51–100 m range, strong in the 0–50 m (${\mathrm{p}}_{\mathrm{s}}$ = 0.61) and 101–150 (${\mathrm{p}}_{\mathrm{s}}$ = 0.69) ranges, and very strong in the 151–200 m range (${\mathrm{p}}_{\mathrm{s}}$ = 0.88). In the SC 2 slope class, a strong correlation was observed in the 0–50 m (${\mathrm{p}}_{\mathrm{s}}$ = 0.77) and 51–100 m (${\mathrm{p}}_{\mathrm{s}}$ = 0.79) ranges, while a very strong correlation was observed in the 101–150 m (${\mathrm{p}}_{\mathrm{s}}$ = 0.81) and 151–200 m (${\mathrm{p}}_{\mathrm{s}}$ = 0.87) ranges. As stated by [90], the closer the value of the coefficient is to one, the stronger the correlation. This allows for the classification of these relationships as either very strong, strong, moderate, weak, or very weak monotonic relationships.

Regarding the productivity and the machine elements of the grapple skidder, a negative correlation was observed for the four machine elements in all skidding distance ranges. In accordance with the findings of [91], one method of establishing relationships between two or more variables under uncertain conditions is the use of correlations. Thus, it was shown that an increase in the time of the machine elements resulted in a decrease in productivity because they tended to move in opposite directions.

For the machine element ET, except for the skid distance range of 101–150 m, where the correlation was significantly weak (0 < |${\mathrm{p}}_{\mathrm{s}}$| ≤ 0.2), the correlations were relatively weak in other ranges (0.2 < |${\mathrm{p}}_{\mathrm{s}}$| ≤ 0.4). For the machine element TL, the correlations were considered weak in the ranges 51–100 m and 151–200 m and moderate in the remaining two ranges (0.4 < |${\mathrm{p}}_{\mathrm{s}}$| ≤ 0.6), as shown in Figure 1.

The machine element BL showed a weak correlation in the skid distance range of 151–200 m, while in the other ranges, it was interpreted as a moderate correlation. For machine element BU, the correlation was found to be weak (0.2 < |${\mathrm{p}}_{\mathrm{s}}$| ≤ 0.4) in the 51–100 m range and very weak (0 < |${\mathrm{p}}_{\mathrm{s}}$| ≤ 0.2) in all other ranges.

In the slope class SC 2, productivity showed a negative correlation with the machine elements of the grapple skidder over all ranges of skid distances. Consequently, the dependence measure showed a similar pattern to that observed in the SC 1 slope class. The authors [92] defined correlation as a valuable numerical statistic for measuring the association between two random variables. In light of these findings, it can be postulated that the highest values for machine element times are associated with the lowest productivity values observed for the grapple skidder.

For the machine element ET, the correlation was judged to be very weak except for the skidding distance range of 51–100 m and weak in all other ranges. For the TL machine element, the correlation was considered weak in all skid distance ranges. Although these correlations were found to be statistically weak, the influence of these variables should not be dismissed in order to avoid inconsistency in estimating the productivity of the grapple skidder. According to [93], the magnitude of an effect should always be evaluated.

For the machine element BL, the coefficient value was found to be moderately correlated only in the 51–100 m range. In all other ranges, the correlation was considered weak. The BU machine element also had the lowest coefficients when compared to the SC 1 slope class values, with the correlation classified as weak for the 0–50 m range and very weak for the remaining skid distance ranges (Figure 2).

In the context of planted forests, these productivity results of the grapple skidder can provide support for improving policies for granting and maintaining forest certification. Most of the planted forests in Brazil are already certified, mainly by the FSC certification system [94]. Better and more transparent forest management explains the movement towards forest certification in Brazil, with certificate holders reporting high overall satisfaction with market access [95]. In addition, there is relatively strong evidence that forest certification has a positive impact on productivity, creating an urgent need for companies to promote and increase productivity to ensure long-term sustainability [96].

3.3. Stochastic Cost Analysis

The mean cost per scheduled hour of the grapple skidder was 83.34 USD h⁻¹ ± 2.29 USD h⁻¹, with a minimum value of 72.77 USD h⁻¹ and a maximum of 95.19 USD h⁻¹. The best fit of the data was found to be the normal probability distribution (Figure 3). Accordingly, these values can be used to meet the budgetary objectives set by the forest enterprise or to set a price for the commercialization of the tree-cutting activity. As stated by [97], these estimates can be subsequently refined as incurred costs become available.

The cost per scheduled hour of a forestry machine depends on a number of factors, including the residual value at the end of its useful life, its economic life, the maintenance conditions, the interest rate used to repay the capital, and the engine power. It should be noted that the inclusion of these factors may or may not be taken into account by the various methodologies used for this purpose. Furthermore, reference [98] states that it is difficult to assess the effectiveness of any particular methodology. However, they emphasize the importance of selecting a method that is appropriate to the characteristics of the equipment in question.

The descriptive statistics for the components of the cost per scheduled hour of the grapple skidder are shown in

Table 5. Descriptive statistics for the components of the cost per scheduled hour (USD h⁻¹) of the John Deere 848H grapple skidder.

Descriptive Statistics	Depreciation	Return on Capital	Insurance	Shelter	Property Taxes	Grapple Skidder Transportation	Fuel	Lubricating Oil and Grease	Maintenance and Repair	Pneumatic	Labor with Operator	Overheads
Minimum	6.60	2.93	0.97	0.24	0.48	0.29	19.86	3.97	8.64	2.21	21.32	3.38
Maximum	8.91	3.97	1.31	0.33	0.65	0.39	26.82	5.36	11.68	2.99	28.83	4.56
Mean	7.76	3.45	1.14	0.28	0.57	0.34	23.34	4.67	10.16	2.60	25.08	3.97
Standard deviation	0.47	0.21	0.07	0.02	0.03	0.02	1.43	0.29	0.62	0.16	1.54	0.24
Asymmetry	0.003	0.003	−0.005	−0.001	0.005	−0.003	−0.001	−0.005	0.002	0.007	0.001	−0.003
Kurtosis	2.40	2.40	2.40	2.40	2.40	2.38	2.41	2.40	2.40	2.39	2.41	2.40
Percentiles
5%	6.96	3.10	1.02	0.26	0.51	0.31	20.94	4.19	9.12	2.33	22.50	3.56
15%	7.23	3.22	1.06	0.26	0.53	0.32	21.76	4.35	9.47	2.42	23.38	3.70
25%	7.42	3.30	1.09	0.27	0.54	0.33	22.32	4.46	9.72	2.48	23.98	3.79
35%	7.56	3.37	1.11	0.28	0.55	0.33	22.77	4.55	9.91	2.54	24.46	3.87
45%	7.70	3.42	1.13	0.28	0.56	0.34	23.16	4.63	10.08	2.58	24.88	3.94
55%	7.81	3.48	1.15	0.29	0.57	0.34	23.52	4.70	10.24	2.62	25.27	4.00
65%	7.95	3.53	1.17	0.29	0.58	0.35	23.91	4.78	10.41	2.66	25.69	4.07
75%	8.10	3.60	1.19	0.30	0.59	0.36	24.35	4.87	10.61	2.71	26.17	4.14
85%	8.28	3.69	1.21	0.30	0.61	0.36	24.91	4.98	10.85	2.78	26.78	4.24
95%	8.55	3.81	1.25	0.31	0.63	0.38	25.72	5.14	11.21	2.87	27.66	4.37

. The asymmetry and kurtosis of these components indicate that the distributions of these components have a pattern close to a normal probability distribution, with respective values of approximately 0 and 3. In other words, the data showed that the frequency of the distribution did not deviate from a symmetrical position, resulting in mesocurtic kurtosis. Consequently, the frequency curves were identical to those of the normal probability distribution, as previously reported by [99,100].

The mean cost of operator labor (25.08 USD h⁻¹) was the most representative among the other costs that constituted the cost per scheduled hour of the grapple skidder, followed by the mean cost of fuel (23.34 USD h⁻¹). This can be attributed to the proportion of social security contributions and benefits applied to the remuneration of the operators. These cost components were also identified as the most important by [101,102] in their analysis of timber harvesting carried out by self-propelled forest machines with similar technical characteristics.

From this perspective, social benefits, such as transportation for operators, meals, and health and dental plans, contribute to the total cost of labor. However, it is important to consider these expenses in the broader context of labor costs. On the other hand, the cost of fuel necessitates the implementation of measures to minimize it, with the aim of reducing the cost per scheduled hour. In other words, as stated by [103], fuel consumption is a function of engine power, which is a challenge in the selection of forestry machines. In addition, operator training and education, with a focus on strategies to reduce fuel consumption, can help reduce the expenditures allocated to this cost component.

The correlations between the cost components and the cost per scheduled hour were positive, indicating that an increase in the value of the components resulted in an analogous increase in the cost per scheduled hour. Accordingly, the cost of operator labor and fuel had the highest Spearman correlation coefficients, ${\mathrm{p}}_{\mathrm{s}}$ = 0.66 and ${\mathrm{p}}_{\mathrm{s}}$ = 0.61, respectively. These values indicate a strong correlation, which supports the findings of [104]. In their study, these costs were identified as the most significant among those evaluated in activities involving the grapple skidder.

The correlation between maintenance and repair costs (${\mathrm{p}}_{\mathrm{s}}$ = 0.26) and depreciation (${\mathrm{p}}_{\mathrm{s}}$ = 0.20) was found to be weak, while the correlation between lubricating oil and grease costs and the remaining cost components was found to be very weak (${\mathrm{p}}_{\mathrm{s}}$ = 0.12). Correlations of less than 0.10 were observed for the remaining cost components. Although these cost components showed weak correlations, they should not be dismissed because a positive correlation would result in an underestimation of the cost per scheduled hour as suggested by [105,106].

Regarding the monotonic relationship between the volume of skidded timber and the production cost, the correlations were negative (p-value < 0.05) in both slope classes SC 1 and SC 2 and in all distance ranges, indicating that a decrease in the volume of skidded timber resulted in an increase in the production cost of the grapple skidder. In SC 1, the Spearman correlation coefficients indicated that the correlations were strong in the 0–50 and 51–100 m ranges (${\mathrm{p}}_{\mathrm{s}}$ = 0.69 and ${\mathrm{p}}_{\mathrm{s}}$ = 0.71, respectively) and very strong in the 101–150 (${\mathrm{p}}_{\mathrm{s}}$ = 0.87) and 151–200 m (${\mathrm{p}}_{\mathrm{s}}$ = 0.88) ranges. In SC 2, the bands 0–50 m (${\mathrm{p}}_{\mathrm{s}}$ = 0.71) and 51–100 m (${\mathrm{p}}_{\mathrm{s}}$ = 0.74) were classified as strong, while the ranges 101–150 m (${\mathrm{p}}_{\mathrm{s}}$ = 0.84) and 151–200 m (${\mathrm{p}}_{\mathrm{s}}$ = 0.86) were classified as very strong.

Regarding the production cost of the grapple skidder in slope class SC 1, the data showed a normal probability distribution for all skid distance ranges (Figure 4). The independent variable accumulated frequency, which ranges from 0 to 1, represents the cumulative sum of the frequencies associated with each production cost interval, taking into account the different skid distance ranges [107]. It reflects the total proportion of scenarios where the cost is less than or equal to the upper limit of each interval along the skid distance ranges.

In SC 1, considering all skid distance ranges, the mean production cost was 0.82 USD m⁻³. For example, in the distance range of 0–50 m (red curve with mean production cost of 0.62 USD m⁻³ ± 0.32 USD m⁻³ and BIC = 24,992.03), 60% of the scenarios have production costs below 0.80 USD m⁻³. In the 51–100 m distance range, the mean production cost was found to be 0.67 USD m⁻³ ± 0.14 USD m⁻³ (BIC = 842,950.87). In the 101–150 m distance range, the mean production cost was 0.87 USD m⁻³ ± 0.20 USD m⁻³ (BIC = 34,992.03), while in the 151–200 m distance range, it was 1.11 USD m⁻³ ± 0.27 USD m⁻³ (BIC = 26,144.08).

In SC 2, considering all skid distance ranges, the mean production cost was 1.48 USD m⁻³. Production costs in the SC 2 slope class (Figure 5), with the exception of the 0–50 m skid distance, showed the best fit when analyzed using the Laplace probability distribution (BIC = 377,715.73) and a mean value of 0.95 USD m⁻³ ± 1.72 USD m⁻³. However, the majority of the remaining data sets showed a better fit when analyzed using the normal probability distribution.

The mean production cost for the 51–100 m distance range was 1.40 USD m⁻³ ± 0.63 USD m⁻³ (BIC = 153,220.93). In the skid distance range of 101–150 m, the mean production cost was 1.66 USD m⁻³ ± 0.74 USD m⁻³ (BIC = 180,507.95). In the 151–200 m distance range, the mean production cost was 1.92 USD m⁻³ ± 0.80 USD m⁻³ (BIC = 196,815.35).

It can be seen that in the classes of flat and wavy to strong wavy slopes, the skidding distance was directly proportional to the production cost of the grapple skidder. This indicates that as the distance in question increased, the production cost also increased. In addition, the greater the distance traveled by the grapple skidder, the longer the travel time, which in turn reduces productivity and is likely to increase fuel consumption. This is a significant cost component, accounting for 26.8% of the total cost of skidding timber by grapple skidder in Eucalyptus planted forests.

4. Conclusions

The mean percentage of the total operating cycle time of the grapple skidder in Eucalyptus planted forests spent on empty travel and travel loaded is 73.8%.

The triangular and exponential probability distributions, respectively, show superior performance in adjusting the times of the grapple skidder machine elements during skidding of timber in the wavy to strong wavy and flat slope classes.

Increasing the skidding distance, which is limited to 200 m in Eucalyptus plantations on a flat slope, reduces grapple skidder productivity by an average of 19.1%, while on a wavy to strong wavy slope the average reduction is 15.1%.

The monotonic relationship between the volume of timber skidded and productivity alters the magnitude of grapple skidder productivity in the same direction, showing 50.0% strong positive correlation, 38.0% very strong positive correlation, and 12.0% moderate positive correlation across the four skidding distance ranges and two slope classes.

The mean productivity of the grapple skidder during timber skidding in Eucalyptus planted forests on flat slopes is 42.2% higher than that observed on the wavy to strong wavy slopes.

Spearman’s correlation coefficients of the monotonic relationship between grapple skidder productivity and machine element times show a negative correlation, indicating a decrease in productivity due to an increase in time spent. This correlation is observed to have a weak to moderate influence.

The mean values for operator labor and fuel are the most significant cost components, together accounting for 58.1% of the cost per scheduled hour of the grapple skidder.

The correlation coefficients show a robust positive relationship between operator labor costs and fuel costs with the cost per scheduled hour, indicating that these variables have a significant influence on the overall composition of grapple skidder costs in Eucalyptus planted forests.

The mean production cost of the grapple skidder for skidding timber in Eucalyptus planted forests on flat slopes is 80.0% lower than on wavy to strong wavy slopes.

In the skid distance range up to 100 m, there is a strong negative correlation between the volume of skid timber and the production cost. This correlation is even more pronounced for skid distances up to 200 m.

By identifying the specific characteristics of each probability distribution and the correlation values in different skidding distance ranges and slope classes, the application of the Monte Carlo method allows more accurate modeling for optimizing the use of grapple skidders in Eucalyptus planted forests.

Some limitations of this study should be considered when interpreting the results. There are unobserved factors (i.e., road structures, species composition, silvicultural management, age structure, etc.) that have a potential impact on grapple skidder productivity and costs that could not be observed.

For further study, it is recommended that the Monte Carlo method be applied to other activities inherent to forest harvesting operations that include other planted forest species of commercial interest, such as Pinus sp. Furthermore, the impact of other slope classes, skid distances, and silvicultural aspects, such as average individual tree volume, on productivity and production costs can be considered.