Productivity and Costs of Mechanized Skidding operations at Sao Hill Forest Plantation, Tanzania

Abstract

Due to global advancement of technology in forest operations, utilization of advanced machineries such as grapple skidder (GS) in timber harvesting has been increasing in the last decades. However, in order to understand their contribution in sustainable harvesting operations, it is important to understand their performance under different operating environment. Therefore, this study aimed to quantify productivity and cost of mechanised skidding operations at Sao Hill Forest plantation (SHFP). Six variables; diameter at a breast height (dbh), tree height, skidding distance, slope, costs, and cycle time (determined using detailed continuous time study) were collected in 120 GS observations.GS productivity and costs were estimated using productive machine hour (PMH) and delays inclusion approach. Regression models were developed using a generalized linear model (GLM) approach. GS productivity under PMH was 2.6% higher than the one including delay time, while skidding costs was 2.1% higher in the approach including delays. This study revealed significant variations (p-value <0.05) in productivity and cost on various terrain classes. At 0 m – 50 m distance, an average delays free GS productivity was 85.5 m³/h, with costs amounting to 1.7 USD/m³. On the distance exceeding 150 m, productivity dropped to 20.1 m³/h, and costs increased to 12.7 USD/m³. Likewise, in 0.0% - 10.0% slope range, average delays free GS productivity and costs was 100 m³/h and 1.5 USD/m³ respectively, while at 20.1% - 30.0% slope range, productivity dropped to 32.6 m³/h and costs raised to 3.9 USD/m³. Skidding distance, slope, and volume per trip were robust predictors of the GS productivity and costs, yielding pseudo-R² values of 58.1% and 64.3%, respectively. Therefore, this study developed statistical models useful for predicting GS productivity and costs, however, their applications are recommended to be within the ranges of the variables used to develop the models.

Keywords:

• cycle time
• delay time
• grapple skidder
• terrain variables
• work element

Introduction

Forest sector particularly plantation forests form a significant natural resource that is renewable in both economic and environmental sense (Spinelli et al., 2020). Therefore, proper management through careful and sustainable utilization in connection with climate change are very essential (Hartsch et al. 2021) and are highly encouraged to ensure sustainability to subsequent generations (Visser and Stampfer, 2003). This can be achieved by ensuring all key operations including timber harvesting are efficiently planned and managed (Parajuli et al., 2020).

To ensure proper planning and management of these operations particularly for timber harvesting so as to meet expected timber quality and demand (Jonsson et al., 2023), it is essential to consider various factors which limits their productivity and cost. Among others include; species type, age, stand density, timber size, bunching strategy, learning curve of the crew (Varch et al. 2021), site conditions (slope, soil, obstacles), available machine types, harvesting prescriptions, intended final products (Akay, 1998; Behjou, 2018; Jiroušek et al., 2007; Kluender and Stokes, 1996; Mizaras et al., 2008; Parajuli et al., 2020), training and motivation level to forest workers (Šporčić et al., 2023).

With the advancement of technology (i.e., due to widespread availability of information technologies (IT) (Llorente et al. 2023; Spinelli et al. 2019), increase in timber demand driven by urbanization and population growth, increase in safety level (Obi and Visser 2020), as well as intentions for costs reduction (Varch et al. 2021), logging operations have gradually shifted toward full mechanization (Mizaras et al., 2008). This transition has been facilitated by the introduction of highly efficient machines such as farm tractors, feller bunchers, and skidders (Bavaghar et al., 2010; Borz et al., 2015; Gölci et al., 2018; Orlovský et al., 2020). These machines are known for their significant production capacity and cost-effectiveness when properly planned. This had driven various researchers across the world to see the importance of analysing productivity and costs of these machines based on real-world scenarios to enhance its proper utilization to its fully capacity. Various techniques, including gross time study, detailed time study (Borz et al., 2013) and regression analysis such as Ordinary Least Square (OLS) have been widely applied in global research (Akay et al., 2004a) to estimate time consumption, productivity (Jonsson et al. 2023) and determine the impact of environmental, stand, economic, and human factors on machinery productivity and costs (Orlovský et al., 2020; Proto et al., 2018). For instance, studies by Ghaffariyan (2020a), Kulak et al. (2017), Sabo and Poršinsky (2005), and Stoilov et al. (2021) applied multiple linear regression analysis on predicting GS time consumption, productivity and costs. Under all reviewed studies, the variables; number of logs, slope, skidding distance, bunching distance and volume per turn were observed to be significant (p-value < 0.05). Moreover, study by Borz et al. (2013) applied Stepwise backward regression techniques which normally involved combining all possible independent variables in the model at the beginning, followed by redundancy of insignificant variables at the end in order to bring highly efficient predictive models.

Although OLS have gained popularity in the field of forest management, its applications is limited by a number of assumptions (Liski et al., 2020), consequently, log and back transformation have always been used to account for the assumptions related with linearity and normality. In light of these limitations, a generalized linear model (GLM) has emerged as a valuable extension to the regression model (Zhang et al., 2020). Unlike traditional methods, the GLM can estimate machinery productivity more flexibly without requiring data transformation. It offers a better fit to the original data, handles collinear variables, and facilitates formulation of highly precise models for accurate estimations of the machinery performance in the field of forest operations (Orlovský et al., 2020).

Despite their predictive usefulness, models for estimating productivity and costs in mechanized timber harvesting operations are limited particularly for the studies conducted in the Sub-Saharan Africa, including Tanzania. Furthermore, productivity and cost of timber harvesting machineries are mostly influenced by environmental factors, stand characteristics, and machine-specific attributes which are highly contingent upon the specific locality and the machine itself. Consequently, there is a pressing need for further research to quantify and develop predictive models that can provide most accurately localized information for future predictions.

In Tanzania mechanized timber harvesting operations begin since early 2000s through integration of various harvesting systems such as; short wood, tree length and whole tree system using various machines such as; feller buncher, cable skidder, GS, farm tractor and three wheeled loaders. Various companies, including Mufindi Paper Mill (MPM), Green Resources Company Ltd. (GRL), and Tanzania Wattle Company (TANWAT) which are situated in the southern highlands, own a majority of these machines (Mauya et al., 2011). Despite of their known production capability, accurate machinery productivity and costs estimations still remain challenging due to largeness of SHFP, its diverse terrain and stand attributes. Therefore, this study addresses this knowledge gap by (i) quantify the productivity and costs of GS using the whole tree harvesting system, (ii) develop models that encompass GS time consumption, productivity, and costs at SHFP (iii) predicting GS productivity and costs across different operating terrain conditions at SHFP using the developed models. It is anticipated that the results from this study will serve as a valuable guide (Akay et al., 2004b; Spinelli et al., 2002) to logging managers, contractors, and forest owners (Miyajima et al., 2021) in accurately estimating and efficiently utilizing harvesting machinery in diverse locations (Ackerman et al., 2018; Borz et al., 2015; Gilanipoor et al., 2012).

Materials and methods

Study area description

This study was conducted at SHFP located in the Southern highland of Tanzania, Iringa region, Mufindi District. The plantation lies between 8° 18′ S to 8° 33′ S latitudes and 35° 06′ E to 35° 20′ E longitudes (Figure 1).

Figure 1. Map of SHFP showing the compartment understudy.

The plantation is about 18km from Mafinga town with estimated total area of 135,903 ha. Due to its largeness, currently the plantation is administratively segmented into four blocks/divisions namely: Irundi, Ihefu, Ihalimba and Mgololo (MNRT, 2018). Each division is broken into compartments differentiated by species, planting date and site class. Generally, the plantation is composed of exotic softwood tree species of Pinus patula, Pinus elliottii, Pinus caribeae, and Cupressus lusitanica, together with hard wood species of Eucalyptus maidenii, Eucalyptus saligna and Eucalyptus camaldulensis. Also, the plantation comprises small patches of natural vegetation characterized by the mosaic of open grassland with scattered indigenous tree species of Erythrina abyssinica, Parinari curatellifolia, Apodytes dimidiata and Albizia petersiana.

Topography and Climate

The plantation is on rolling terrain interacting with some low hills and wide flat-bottomed valleys with an altitude ranging from 1,700m to 2,000 m above sea level. The soils of the area are acidic deep soil with a pH range from 4.4 to 5.4. The climate is characterized by unimodal rainfall pattern starting from November to April and a dry season from May to late October. The mean annual rainfall is 1,300 mm with a range of 725 mm to 1,400 mm. The temperatures are fairly cool, with the minimum monthly temperature range between 10 °C to 18 °C while the maximum range is between 23 °C to 28 °C (MNRT, 2018).

The study focused on mechanized timber skidding operations using the whole tree system. Harvesting operations in the study area involves clear-cutting, utilizing both semi-mechanized and mechanized approaches. Mechanized operations employ machinery like feller bunchers, farm tractors with wire rope winches, cable and grapple skidders. Semi-mechanized operations primarily rely on chainsaws and farm tractors. The main harvesting systems employed are cut-to-length, tree-length, and whole tree systems.

Sampling procedures

The study employed two phase sampling. On the first phase, forest compartment for the study was obtained with the aid of SHFP management plan in order to avoid destruction of immature forest stands. Moreover, since our study focused on ground skidding operation, terrain variables were considered as one of the factors that influence the performance of the GS. Therefore, forest compartment named ILASA 3/9 with heterogeneous terrain was selected with the aid of digital elevation model (DEM) with 30m x 30m resolution downloaded from https://earthexplorer.usgs.gov/. A total of 30 square plots of 15 m x 15 m were randomly distributed throughout the compartment to capture all ranges of terrain (Figure 2). On second phase, 120 GS work cycle with respective work elements was determined using the formula by Murphy, (2005) as indicated in equation 1(1) $$\mathrm{n}={t^2\mathrm{*}\mathit{\text{Var}}\left(\mathit{\text{WCT}}\right)/[E\mathrm{*}\overline{\mathit{\text{WCT}}}/100]}^2$$(1) . Whereby, 10 pilot time study was conducted, giving an average cycle time (i.e., mean WCT) of 6.873 minutes, Variance cycle time (i.e., Var WCT) of 0.071 minutes and 0.95 precision level (E) was used.

Figure 2. Sampling procedure used to select forest compartment understudy.

(1) $$\mathrm{n}={t^2\mathrm{*}\mathit{\text{Var}}\left(\mathit{\text{WCT}}\right)/[E\mathrm{*}\overline{\mathit{\text{WCT}}}/100]}^2$$(1)

Where; n = number of work cycles to be studied, t = Student’s t-value, Var (WCT) = Variance of the work cycle time, E = Level of precision desired, $\overline{\mathrm{\text{WCT}}}$= Mean work cycle time (minutes).

Machinery description

The machinery used under this study was CAT 525 grapple skidder, purchased by MPM in 2018. It is American manufactured caterpillar brand with the following characteristics available in www.ritchiespecs.com as presented in Table 1.

Table 1. GS descriptions

Machinery specifications	Description
Model	CAT 525
Configuration	Rubber-tired
Overall length	6.00 m
Overall width	3.13 m
Ground clearance	0.53 m
Wheel base	3.50 m
Engine Model	CAT 3304 DIT Diesel
Gross power	175 hp
Operating weight	15200 kg
Brakes service type	Hydraulic actuated, oil disc
Maximum Drawbar pull	19731.3 kg
Maximum forward speed	16.9 mph
Maximum reverse speed	12 mph
Estimated operating weight	13558.3 kg
Grapple bunching capacity	1.161m^2 (1065.9 kg)
Maximum operating distance	500 m
Maximum operating slope	30 %

CAT = Caterpillar, hp = Horse power, kg = Kilogram, m = meter, mph = Miles per hour.

Data collection

GS cycle time (minutes), individual tree variables (dbh and height), terrain variables (skidding distance and slope) and costs variables were collected. The details for each collected variable are described below.

GS cycle time

GS productive time (PMH) which refers to the time spent by the skidder to perform a given task at the workplace (Magagnotti et al., 2013) was recorded through detailed continuous time study technique using stopwatch. The GS cycle time was segmented into five work elements namely: travel empty (TE), positioning (PS), grappling (GL), travel loaded (TL) and unloading (UNL). Also, necessary delays (inevitable interruption due to the nature of the work and the environment), as well as unnecessary delays (the wastage of time which can be eliminated by improving supervision and training to workers) (Mauya 2022), were recorded on each work element of the GS once they occurred. Time recording started at the beginning of a work cycle and the stopwatch was paused at the end of the cycle. Elapsed time was read directly from the stopwatch and recorded on the prepared field forms.

Tree variables

Before the beginning of skidding operations, all standing trees within 15m x 15m square plots were labelled by numbers. Thereafter, dbh and tree height for each tree were measured using calliper and vertex hypsometer respectively. Since skidding operation was performed using whole tree harvesting system (i.e., with branches and tops intact to the main trunk), the volume of each tree was estimated using the allometric single tree equation (Equation 2(2) $$\text{Tree volume}=\text{exp}(-9.04925+\text{1}.14781\times \text{ln}\left(h\text{eight}\right)+\text{1}.5496\times \text{ln}\left(\text{dbh}\right)$$(2) ) by Malimbwi et al., (2016). The model applied is compatible with Pinus patula grown in the study area. (2) $$\text{Tree volume}=\text{exp}(-9.04925+\text{1}.14781\times \text{ln}\left(h\text{eight}\right)+\text{1}.5496\times \text{ln}\left(\text{dbh}\right)$$(2)

Terrain variables

Skidding distance (m) and slope (%) were the two-terrain variables recorded during skidding operation. Skidding distance from the stump site to the landing was measured using measuring tape, while slope was measured using vertex hypsometer.

GS costs variables

GS operational (i.e., fixed and variable) costs information was collected by interviewing logging supervisor and procurements personnel from Mufindi Paper Mill. Simply, fixed (standing or capital) costs are the ones that need to be recovered by machine irrespective of the amount of work a machine does or the revenue it earns and are associated only with owning the machine. Normally includes: depreciation, interest, insurance, and taxes (Ackerman et al., 2014). On the other hand, variable (running) costs are the ones incurred when the machine is working, travelling empty, or when the engine is running (Bjorheden and Thompson, 1995). They normally vary directly with the level of output produced by the machine. It includes costs for labour, fuel, lubricants, machine maintenance and tires. (Ackerman et al., 2014; Nwokoye and Ilechukwu, 2018).

Data analysis

GS productivity (m³/h) was computed using two approaches; with inclusion and exclusion of necessary delays i.e., equations 3(3) $${\mathrm{P}}_{\mathrm{D}}=\frac{\mathrm{\text{Tvol}}\mathrm{ }\left({\mathrm{m}}^3\right)\text{x F x}60}{\mathrm{T}}$$(3) and 5(5) $$\mathrm{P}=\frac{\mathrm{\text{Tvol}}\mathrm{ }\left({\mathrm{m}}^3\right)\mathrm{x}\left(60\right)}{\mathrm{\text{PMH}}}\mathrm{ }$$(5) respectively. (3) $${\mathrm{P}}_{\mathrm{D}}=\frac{\mathrm{\text{Tvol}}\mathrm{ }\left({\mathrm{m}}^3\right)\text{x F x}60}{\mathrm{T}}$$(3)

Where: P_D is the GS productivity (m³/h) per work cycle (including necessary delay), Tvol is skidded tree volume (m³) per work cycle, 60 is the number of minutes per a given workplace hour and T is total GS productive time (minutes) per work cycle, F is the proportion of productive time. (4) $$\mathrm{F}=\frac{100‐\mathrm{D}}{100}$$(4)

Where D is Delay time expressed as a percentage of workplace time in minutes.

On the other hand, GS productivity without delay time P_(PMH) was computed using the formula adapted in the studies by Mauya (2022), and Miyajima et al., (2021a) Equation 5(5) $$\mathrm{P}=\frac{\mathrm{\text{Tvol}}\mathrm{ }\left({\mathrm{m}}^3\right)\mathrm{x}\left(60\right)}{\mathrm{\text{PMH}}}\mathrm{ }$$(5) . (5) $$\mathrm{P}=\frac{\mathrm{\text{Tvol}}\mathrm{ }\left({\mathrm{m}}^3\right)\mathrm{x}\left(60\right)}{\mathrm{\text{PMH}}}\mathrm{ }$$(5)

Where: P is the GS productivity (m³/h), 60 is the conversion factor for converting minutes to hours and PMH is the productive machine hour in minutes.

Lastly, paired t-test was performed to test if there is a significant difference between GS productivity with delay and without delay time.

All GS cost components (Equations 6(6) $$\mathrm{\text{Depreciation}}\left(\mathrm{D}\right)=\frac{\text{Delivered cost x}0.90}{\text{life time}\left(\mathrm{\text{hours}}\right)}$$(6) to 16(16) $$\text{Unit skidding costs}(\mathrm{\text{USD}}/{\mathrm{m}}^3)=\frac{\text{Hourly skidding costs}(\mathrm{\text{USD}}/\mathrm{h})}{\text{Production rate}({\mathrm{m}}^3/\mathrm{h})}$$(16) ) were computed based on the formulas described by Sessions, (2007). The machine (i.e., GS) under study was CAT brand, with registration number T.809 DNB, having 175 horse power and it was purchased in the year 2018.The machinery delivered costs was 304,117.4 USD and has saved over a period of six years up to the present study. Furthermore, the machinery annual tax was amounting to 24,472.3 USD/year, with zero interest rate since the machine was purchased cash. The average working time was 3506 h/year, with the labour (operator’s) costs amounting to 0.56 USD/h.

Fixed costs

1. Depreciation (6) $$\mathrm{\text{Depreciation}}\left(\mathrm{D}\right)=\frac{\text{Delivered cost x}0.90}{\text{life time}\left(\mathrm{\text{hours}}\right)}$$(6)

2. Insurance (7) $$\mathrm{\text{Insurance}}\left(\mathrm{I}\right)=\frac{\text{Delivered costs x}0.60\mathrm{ }\mathrm{x}0.03}{\text{Average hours per year}}$$(7)

3. Interest (i) (8) $$\mathrm{\text{Interest}}\mathrm{ }\left(i\right)=\frac{\text{Delivered cost x}0.60\text{x interest rate}\left(\mathrm{\%}\right)}{\text{Average hours per year}}$$(8)

4. Taxes (Tc) (9) $$\mathrm{\text{Taxes}}\left(\mathrm{\text{Tc}}\right)=\frac{\text{Anual tax amount}}{\text{Average hours per year}}$$(9)

Variable costs

1. Labor costs (Lc) (10) $$\text{Labor cost}\left(\mathrm{\text{Lc}}\right)=\frac{\text{Labour salary per period x}(1+\mathrm{f})}{\text{Machine hours in period}}$$(10) Where; f = Fringe benefits expressed as % of direct labor costs (i.e., 0.5% in this case).

2. Fuel cost (Fc) (11) $$\text{Fuel costs}\left(\mathrm{\text{Fc}}\right)=\mathrm{\text{GHP}}\mathrm{ }x\text{X}x\mathrm{ }\mathrm{\text{CL}}$$(11) Where: GHP = Gross engine horsepower (i.e., 175 HP in this case), CL = Fuel price per liter (i.e., 1.236 USD/liter in this case), X = 0.12 for diesel and 0.175 for the gasoline engine.

3. Lubricants costs (Oil and grease) (12) $$\text{Lubricant costs}\left(\mathrm{L}\right)=\frac{\mathrm{\text{GHP}}\mathrm{ }x\text{X}x\mathrm{ }3.4}{100}$$(12) Where; GHP = Gross engine horsepower (i.e., 175 HP in this case), and X = 0.20 for tractors, skidders, front end loaders and trucks, X = 0.30 for feller bunchers and knuckle boom loaders, X = 0.50 for processors, harvesters and forwarders.

4. Maintenance costs (Mc) (13) $$\text{Maintenance costs}\left(\mathrm{\text{Mc}}\right)=\frac{\text{Delivered costs}}{\text{Life time}\left(\mathrm{\text{hours}}\right)}$$(13)

5. Tires costs (Ty) (14) $$\text{Tires costs}\left(\mathrm{\text{Ty}}\right)=0.0006\text{x CST}$$(14) Where; CST = Price of a set of the replaced tires (i.e., 2,606.7 USD in this case).

Hourly skidding costs (USD/h)

(15) $$\text{Hourly skidding costs}\left(\text{USD}/\mathrm{h}\right)=\text{D}+\text{I}+I+\text{Tc}+\text{Lc}+\text{Fc}+\text{L}+\text{Mc}+\text{Ty}$$(15) (16) $$\text{Unit skidding costs}(\mathrm{\text{USD}}/{\mathrm{m}}^3)=\frac{\text{Hourly skidding costs}(\mathrm{\text{USD}}/\mathrm{h})}{\text{Production rate}({\mathrm{m}}^3/\mathrm{h})}$$(16)

Statistical analysis

Model development

Previous studies (e.g. Conrad et al., 2013, Long, 2003, Wang et al., 2004, 2005) were used as the baseline for developing GS time consumption, productivity and costs models for this study. The normality of the data set was tested using Shapiro wilk test whereby variables with a p-value less than 0.05 were considered to be not normally distributed (Yap and Sim, 2011; Korkmaz et al., 2014). Modelling was performed using a GLM approach (Equation 17(17) $$g\left(\mathrm{\mu }i\right)=\mathit{\text{ni}}=\mathrm{\alpha }+{\mathrm{\beta }}_1{\mathrm{x}}_{i1}+{\mathrm{\beta }}_2{\mathrm{x}}_{i\mathrm{ }2}+\dots {\mathrm{\beta }}_k{\mathrm{x}}_{\mathit{\text{ik}}}+\text{ᶓ}$$(17) ) in R-statistical software since the GLMs provide greater flexibility in analysing data even those related to non-normal distributions (Ravindra et al., 2019). Due to the nature of the data being continuous, the Gaussian family was employed to ensure error terms are equally distributed throughout the given data set. Predictor variables for each model were further examined using variance inflation factor (VIF) to test for multicollinearity. Variables with VIF>5 were excluded from the model, indicating a sign of multicollinearity (Shabani et al., 2021). (17) $$g\left(\mathrm{\mu }i\right)=\mathit{\text{ni}}=\mathrm{\alpha }+{\mathrm{\beta }}_1{\mathrm{x}}_{i1}+{\mathrm{\beta }}_2{\mathrm{x}}_{i\mathrm{ }2}+\dots {\mathrm{\beta }}_k{\mathrm{x}}_{\mathit{\text{ik}}}+\text{ᶓ}$$(17)

Where: g(µ_i) = n_i is a smooth and invertible linearizing link function g (·), which transforms the expectation of the response variable, μi = E(Yi), to the linear predictor, α is the model intercept while “β₁ x_{i1 +} β₂ x _{i2 +}….β_k x_{i k”} are independent variables which represent the change in the dependent variable associated with a one-unit change in each independent variable (Lindsey, 1998; Nelder and Wedderburn, 1972).

Finally, pseudo-R-squares (Equation 18(18) $$\text{Pseudo R}‐\mathrm{\text{squared}}=1-\left(\frac{\text{Residual deviance}}{\text{Null deviance}}\right)\mathrm{*}100\mathrm{ }$$(18) ) and residual plots were used to assess the goodness of fit for the developed model. (18) $$\text{Pseudo R}‐\mathrm{\text{squared}}=1-\left(\frac{\text{Residual deviance}}{\text{Null deviance}}\right)\mathrm{*}100\mathrm{ }$$(18)

Model validation

For the precision and accuracy of the developed model to be known, model validation is inevitable (James et al., 2013; Jimmy et al., 2013; Mauya, 2022). Based on the available data set (i.e.,120 GS observations), were segmented into ten-folds randomly. One subset was held out for checking the model performance while the model is trained on the remaining subsets (James et al., 2013). Then residual root mean square error (RMSEr) as an indicator for assessing model quality was determined using its predicted values (Equations 19(19) $$\mathit{\text{RMSE}}=\sqrt{{\sum }_{i=1}^n\frac{{(y}_i-{\mathrm{y}̂}_{i\mathrm{ }}{)}^2}{n}}$$(19) and 20). (19) $$\mathit{\text{RMSE}}=\sqrt{{\sum }_{i=1}^n\frac{{(y}_i-{\mathrm{y}̂}_{i\mathrm{ }}{)}^2}{n}}$$(19) (20) $${\mathrm{\text{RMSE}}}_{\mathrm{r}}=\frac{\mathrm{\text{RMSE}}}{\mathrm{y}̅}\mathrm{x}100$$(20)

Where: ${\sum }_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{\text{is}}$ the sum of all observations from i = 1 to i = n, yi and ŷ denote observed and predicted variables for both time consumption, productivity and unit skidding costs in a given i^th observation respectively.

y̅ denotes the mean for observed time consumption, productivity and units skidding costs for all GS cycles, n is the total number of observations in the dataset.

Results

GS summary statistics

A total of 120 GS work cycles were observed for the entire study. The effective GS workplace time (PMH) ranged from 0.593 minutes at a minimum skidding distance of 14.2 m and slope of 0.4% to 11.878 minutes at a maximum skidding distance of 220 m and slope of 26.0%. The average GS time was 4.519 minutes under PMH approach and 4.625 minutes when including necessary delays. The most GS time-consuming work element was travel loaded (TL) which consumed a substantial part of total skidding time for about 47.9%, while the least work element was unloading (Figure 3). The summary for all variables is presented in Table 2.

Figure 3. GS time distribution for each work element at a mean distance of 59.2 m and slope of 13.5%.

Table 2. Skidding summary statistics

	Variable	Statistical parameter
		Minimum	Maximum	Mean	Standard deviation
Dependent variables	Travel empty (min)	0.180	4.817	1.308	1.003
	Positioning (min)	0.003	2.150	0.576	0.542
	Grappling (min)	0.010	1.483	0.387	0.318
	Travel loaded (min)	0.083	6.117	2.165	1.314
	Unloading (min)	0.033	0.250	0.084	0.043
	Necessary delay (min)	0.000	1.700	0.106	0.256
	Unnecessary delay (min)	0.000	1.633	0.014	0.149
	Total productive time (min)	0.593	11.878	4.519	2.498
	Total time including delays (min)	0.593	11.878	4.625	2.516
Independent variables	Tree dbh (cm)	12.000	69.000	32.497	11.581
	Tree height (m)	11.600	33.600	23.555	4.590
	Tree volume (m^3)	0.123	3.807	1.086	0.716
	Volume per trip (m^3)	0.259	8.352	3.890	1.770
	Slope (%)	0.400	26.600	13.522	7.996
	Distance (m)	14.200	220.000	59.218	36.384

GS productivity and costs

GS productivity for both productive machine hour (PMH) and when including necessary delays, were 67.5 m³/h and 64.1 m³/h respectively. Individual cost items and unit skidding costs for both approaches are presented in Table 3. Furthermore, the performed paired t-test indicated a significant difference between productivity and costs while using the effective time (PMH) and when including necessary delays (p-value < 0.05).

Table 3. GS operational costs

	Value under PMH	Values when including NED
ITEM	Average hourly costs (USD/h)	Average hourly costs (USD/h)
Depreciation	11.73	11.98
Insurance	1.41	1.44
Interest	0.00	0.00
Labor	0.51	0.52
Fuel	70.15	71.67
Lubricants	6.78	6.93
Maintenance	5.86	5.99
Tires	1.41	1.44
Unit skidding cots (USD/m^3)	3.36	3.50

PHM = Productive machine hour, NED = necessary delay time

GS time consumption, productivity and costs models

Time consumption models

The explanatory (independent) and response (dependent) variables vary between models and their number ranged from one to three. The parameter estimates for all models except unloading were significantly different from zero (p <0.05) and the VIF values were <5 indicating an acceptable level of multicollinearity. For the case of time consumption, model for predicting machine productive hours (PMH) and individual work elements i.e., travel empty, positioning, grappling, travel loaded and unloading were developed (Table 4). Pseudo-R-square value which explains the goodness of fit of the model varied among productive work elements. The highest pseudo-R-squared was observed in the total PMH model (81.5%), while the least pseudo-R-squared was observed in unloading (UNL) work element (9.4%). The performance of each model was further assessed through residual plots (Figure 4).

Figure 4. Residual plots for all GS work elements and total productive time.

Table 4. GS time consumption model for individual work elements and total productive time.

Work element	Time consumption model	Pseudo R-squared (%)	RMSEr (%)
Travel empty	−0.069 + 0.012 Slope + 0.02SkD	67.1	37.1
Positioning	−0.116 + 0.051 Slope	57.0	60.4
Grappling	−0.184 + 0.147 Volume per trip	66.7	46.5
Travel loaded	0.320 + 0.028 Slope+ 0.025 SkD	58.0	39.4
Unloading	0.076 + 0.002 Trees per load	9.4	49.3
Total productive time	0.761 + 0.088 Slope + 0.045 SkD - 0.003 Av.trip dbh	81.5	24.7

SkD = Skidding distance, Av.trp dbh = Average dbh per trip

Cross-validation results showed that the RMSEr value for the total productive time model was relatively smaller as compared to individual work elements. However, the scatter plots for positioning and unloading time had shown a behaviour of over-prediction since some of the scatter points deviate far from the average value (Figure 5).

Figure 5. Relationship between observed and predicted GS time on various work elements

Productivity and costs models

For both GS productivity and unit skidding costs, the variables skidding distance (m) and slope (%) appeared to be good predictors for estimating machine productivity and costs. The pseudo-R-squared for productivity and costs models (Equations 21(21) $$\mathrm{\text{Productivity}}({\mathrm{m}}^3/\mathrm{h})=72.651-2.309\mathrm{\text{Slope}}-0.332\mathrm{S}\mathrm{\text{kD}}+10.890\mathrm{A}\mathrm{v}.\text{trip volume}$$(21) and 22(22) $$\text{Unit skidding costs}(\mathrm{\text{USD}}/{\mathrm{m}}^3)=-0.009-0.027\mathrm{\text{Slope}}+0.065\text{SkD}$$(22) ) were 58.1% and 64.3% respectively indicating a good fit for both models. The residual plots for both models are normally distributed, with fewer highly deviated scatter points, which are normally caused by the few observations with higher values (Figure 6). (21) $$\mathrm{\text{Productivity}}({\mathrm{m}}^3/\mathrm{h})=72.651-2.309\mathrm{\text{Slope}}-0.332\mathrm{S}\mathrm{\text{kD}}+10.890\mathrm{A}\mathrm{v}.\text{trip volume}$$(21) (22) $$\text{Unit skidding costs}(\mathrm{\text{USD}}/{\mathrm{m}}^3)=-0.009-0.027\mathrm{\text{Slope}}+0.065\text{SkD}$$(22)

Figure 6. Residual and scatter plots for GS productivity and costs models

The cross-validation results indicated that RMSEr value for the GS productivity and unit skidding cots models were 49.3% and 50.4% respectively. Furthermore, RMSEr were observed to vary across terrain classes. Both productivity and costs models showed lowest RMSEr value at a distance and slope class of 0.00-50.0 m and 0.00 – 10.0% respectively, while the highest RMSEr value was observed at the highest distance and slope class (Table 5). This implies that skidding distance and slope were the main predictors for the GS productivity and costs (Figure 7a and 7b).

Figure 7. a. Error bars show the performance of productivity and cost models on various distances.

b. Error bars show the performance of productivity and cost models on various slopes.

Table 5. Performance of productivity and cost models on various distances and slope classes.

	Skidder Productivity (m^3/h)		Unit skidding costs (USD/m^3)
Distance class (m)	RMSE	RMSEr (%)	RMSE	RMSEr (%)
0.00-50.0	12.1	18.8	0.6	17.1
50.1-100.0	15.4	24.1	1.5	43.4
100.1-150.0	18.0	28.1	2.3	65.9
Slope class (%)	RMSE	RMSEr (%)	RMSE	RMSEr (%)
0.00 - 10.0	40.2	62.7	0.6	18.6
10.1 - 20.0	40.4	62.9	1.5	43.4
20.1 - 30.0	40.7	63.5	2.4	68.6

Discussion

This study aimed to quantify and develop predictive models for the GS time consumption, productivity, and cost using the whole tree harvesting system. The objective was to establish baseline information to accurately estimate and plan mechanized timber harvesting in Tanzania plantation forests and other geographical areas which bare similar environmental and stand characteristics. Regression and correlation analyses were employed to evaluate the influence of machine, environmental, and stand characteristics on GS performance.

The results revealed that the average PMH for the GS was 4.519 and 4.625 minutes, inclusive of necessary delays. The study also identified travel loaded (TL) as the most time-consuming work element, with an average duration of 2.165 minutes, followed by travel empty (TE) i.e., 1.308 minutes. Conversely, unloading time was found to be the least time-consuming element, lasting for only 0.084 minutes. Previous studies by Orlovský et al. (2020) and Borz et al. (2015) also highlighted travel loaded as the mostly time-consuming work element for the GS. The GS time in our study varies with the one reported by Mauya et al. (2011) who conducted similar study in the same forest plantation using John Deere 648G (i.e., 4.8 minutes). By comparing similar aspects which was assessed in both studies, such variations are probably due to the type of machine used, terrain condition and operator experience. This can be justified by higher delay time in the study by Mauya et al. (2011) which took 20.4% of the entire GS productive time mainly due to presence of obstacles in the skidding trails and poor grappling of logs, compared with the one reported in our study (which took only 2.3% of the GS productive time).

The average GS productivity and costs under PHM approach (i.e., 67.5 m³/h and 3.363 USD/m³ respectively) was 2.6% higher than the one including delays, while costs under PMH was 2.1% lower than the one including delays. This was further justified by the results from the paired t-test which indicated that there was statistically significant differences (p-value < 0.05) of productivity and costs between the two approaches. This implies that, through ensuring effective supervision, we can minimise operational delays, improving GS productivity by 2.6% and reduce the skidding operational costs by 2.1%. The GS productivity and costs values under this study fall within the range reported in other studies. For example, Ghaffariyan (2020a) and Klepac and Mitchell (2016) reviewed skidder productivity in coniferous forest plantations and reported a productivity range of 9.3 m³/h to 78 m³/hr. Dodson et al. (2006) reported an average GS productivity of 88.7 m³/h when using the whole tree system in Central Oregon for Juniperus occidentalis forest stands. Adebayo et al. (2007) reported an average productivity of 58.3 m³/h for Cat D-518 GS in northern Idaho, with an average skidding distance of 130 m and slope ranging from 3% to 34%. These studies revealed that average skidded volume per trip, extraction distance, and slope are major factors influencing GS productivity and costs. For instance, Ghaffariyan (2020b) reported lower GS productivity on steep slopes (25% to 40%), ranging from 20 m³/h to 48 m³/h. Additionally, Dodson et al. (2006) reported a unit skidding cost of 1.52 USD/m³ using the conventional timber harvesting system. Furthermore, Akay et al. (2004) reported an average productivity and costs of 19.2 m³/h and 15.35 USD/m³ using CAT 525 GS in Turkish forests when harvesting cedars, pines and firs tree species. The average load volume per turn was 0.68 m³ and the mean slope was 31%. Kulak et al. (2017) reported the GS productivity range of 8 m³/h to 14 m³/h in Pinus sylvestris, southern Poland using John Deere 548G-III in a skidding distance ranged between 124 m and 246 m, and the load volume per turn ranged from 1.86 m³ to 2.95 m³. Moreover, study by Ngulube, (2012) reported an average GS productivity of 21.5 m³/h when using CAT 525C grapple skidder and whole tree harvesting system. The reported terrain slope under this study varied between 5% and 25%, with an average of 11.5%, while the average load volume and skidding distance was 1.45 m³ and 134 m respectively. However, most of these studies were conducted in the geographical regions outside Tanzania and therefore, comparison was aimed on getting the global overview trend of GS productivity and costs. But to our understanding this study is one among the few studies conducted in the plantation forests of Tanzania and the first study to explore the productivity and costs of GS when skidding using whole tree harvesting system, and predicting GS performance using GLM approach.

To enhance applications of GS in the areas with similar terrain and stand characteristics, statistical models for predicting GS time consumption, productivity and costs were developed. The results showed that, all selected predictor variables for the GS time consumption, productivity, and costs predictions were significant (p-value < 0.05), except for unloading (UNL).

The performance of the time consumption models for the different GS work elements was assessed using pseudo-R-square and RMSEr value. The highest pseudo-R-square value (i.e., 81.5%) was observed in the GS total productive time, while the lowest pseudo-R-square value (i.e., 9.4%) was observed in unloading (UNL). The residual plot for unloading showed a funnel pattern, likely due to fewer observations having higher deviation from the average value. A higher pseudo-R-square (i.e., 81.5%) indicated a better fit of the predictor variables in the model, meaning that selected variables; skidding distance, slope, and average tree dbh per trip were the good predictors of the GS skidding time. This is in line with the previous studies by Hiesl (2013), Kopseak et al. (2021), Ngulube et al. (2014), and Vitorelo (2011) which also highlighted skidding distance, slope, and load volume as the main predictors of GS cycle time (p-value < 0.05). For the productivity and costs models, the pseudo-R-square values were 58.1% and 64.3%, respectively, falling within the range reported in other studies e.g., (Visser and Stampfer 2003) reported a GS productivity model with an R-square value of 62%, significantly impacted by skidding distance and volume per load (p-value < 0.05). To further explore the performance of the models, K-fold cross validation was used to test the reliability of the models. The results revealed that prediction accuracy of the model decreased as the terrain variables (i.e., slope and distance) increase (Table 5). Higher GS productivity and costs prediction accuracy (i.e., RMSEr = 18.8% and RMSEr = 17.1%) were observed in a distance and slope class of 0.0 – 50 m and 0.0 – 10.0% respectively, while lower prediction accuracy was found at a higher distance and slope class (p-value < 0.05). Moreover, slight variation in RMSEr value was observed in all slope classes as compared to higher variation among distance classes. That is due to lower variation (Çalişkan, 2019) in the range of the measured filed slope variable (i.e., 13 ± 7.9%) compared to field measured distance (i.e., 59 ± 36 m). Moreover, model’s prediction accuracy in this study may be influenced by other factors which were beyond the scope of this study, including modelling techniques. For instance, study by Munis et al. (2023) showed that other machine learning techniques such as Random Forest has higher prediction accuracy (i.e., RMSE = 18.37 ± 13.10) compared to Linear regression model (i.e., RMSE = 28.82 ± 12.88) when modelling forwarder performance in Brazilian Pinus and Eucalyptus planted forests.

Generally, the findings of this study generated a baseline information in operationalisation of the GS in Tanzania forest plantation, of which such information (i.e., skidding using whole tree harvesting system) was not present before. Moreover, further studies on heterogeneous landscape with more terrain variation are recommended to explore its influence on GS performance. Also, other modelling techniques such as machine learning are encouraged to assess their prediction accuracy on GS performance. Lastly, this study was conducted during the dry season. Hence, these models will be appropriate once machine work on the area bearing similar spatial and temporal characteristics. Further studies of the GS performance on wet condition are also recommended.

Conclusion

The study examined the productivity and costs associated with GS in two different approaches: using PMH and including delay time. The developed models for predicting GS time consumption, productivity, and costs using whole tree conventional timber harvesting system, indicated that variables; distance, slope, average tree dbh, and volume per trip were the significant predictors of the GS performance. These findings provide valuable guidelines to planners for an effective resources and facilities allocation, enabling the fulfilment of timber demand within targeted timelines and maximizing profit. Moreover, the study's findings have broader implications to forest plantations with similar environmental and stand characteristics. The insights gained can contribute to better planning of future harvesting operations, ensuring improved efficiency and productivity of the GS. Further studies are encouraged to explore potentials of other modelling techniques on GS time consumption, productivity and costs predictions.

Author Contributions

Conceptualization: Gilberth Prosper Temba, Ernest William Mauya, George Ansigar Migunga.

Data analysis: Gilberth Prosper Temba.

Methodology: Gilberth Prosper Temba, Ernest William Mauya.

Supervision: Ernest William Mauya, George Ansigar Migunga.

Writing – original draft: Gilberth Prosper Temba.

Writing – review and editing: Gilberth Prosper Temba, Ernest William Mauya, George Ansigar Migunga.

Acknowledgments

We acknowledge the project titled “Higher Education for Economic Transformation (HEET)”, under Sokoine University of Agriculture and Tanzania Forest Fund (TaFF), for providing financial support during data collection.

Additional information

Funding

Sokoine University of Agriculture 10.13039/501100023918 and Tanzania Forest Fund.