Steep Slope Harvest System Models for Small to Large Trees

Background: Tethered cut-to-length and cable yarding systems with tethered falling equipment are increasingly used to harvest trees from slopes exceeding 30–60% more safely and at reduced financial cost than less mechanized harvest systems. Existing studies of harvest equipment typically isolate one or two pieces of equipment in a harvest system and often occur on sites with slopes below 50% and trees less than 60 cm in diameter. Methods: We analyzed machine capabilities and productivity regressions to extrapolate existing models to steep slope harvesting of trees up to 115 cm diameter. The resulting individual machine models are integrated into models of cut-to-length and long-log harvest system productivity. We estimated the financial operating costs of the harvest systems considered from equipment pricing and operator wages. Results: Analysis of even-age Douglas-fir (Pseudotsuga menziesii) and western hemlock (Tsuga heterophylla) rotations suggests eight-wheel forwarder productivity, swing yarder productivity, and mechanization of manual chainsaw labor with tethered harvesters as primary controls on harvest costs. Conclusions: The proposed model enables predictions across a greater range of slopes and tree sizes than those previously modeled, creating a foundation for future research into the cost and productivity of steep slope harvesting systems.


Introduction
Models of timber harvest systems are used to support silvicultural decision-making. Silvicultural optimization systems (e.g., [1,2]) rely on harvest system models, along with other models, to guide land management planning. An inaccurate harvest model can lead to a poor plan and, if limits on harvest equipment capabilities are incorrectly represented, may result in infeasible plans. To implement a land management plan, forest harvest contractors need to determine if equipment in the harvest systems they operate is suitable for a given timber sale. Contractors also need to accurately predict stump to mill costs to bid sales competitively without incurring a financial loss. Harvest system models thus facilitate both land management and harvest contracting by describing the abilities, productivity, operating costs, and interactions of the individual machines in each harvest system of interest. Given an understanding of the trees offered in a timber sale and the logs that would be produced from those trees, a harvest cost model predicts each machine's utilization and the total financial cost for the harvest. Comparing cost estimates between models of different harvest systems helps identify the most appropriate harvest method based on ecological and financial considerations.
On harvest units with steep slopes, cable yarding systems have been used to remove felled trees, or logs cut from felled trees, since the late 1800s [3]. When a slope is too steep for mechanized harvesting equipment, trees are manually felled using chainsaws [4]. One or more logs are often bucked from these trees before yarding, also using chainsaws [5]. Because mechanization reduces risk to human operators and increases productivity compared to manual felling and bucking [6], commercial winch-assist systems enabled use of Figure 1. Stump to mill flow of wood in three steep slope harvest systems. In system (a), a tethered harvester-forwarder pair provides cut-to-length harvesting suitable for thinning and continuous cover forestry, variable retention harvests, or clearfelling. Systems based on yarders (b,c) can be used for thinning but are most efficient in clearfelling [5].
Each piece of heavy equipment in a harvest system has size or weight restrictions on the wood it can handle (Section 2.1). Trees that are too large will not fit in felling or processing heads for cutting or may weigh too much for equipment booms to lift or drag. Logs too heavy for a forwarder's crane cannot be lifted into its bunk. Similarly, skyline deflection determines maximum payload weights in yarding. However, growing large trees on slopes of 50-100% is a frequent silvicultural objective within western North America's temperate rainforests for reasons such as conserving biodiversity, restoring wildlife habitat, carbon sequestration, and raising the financial value of timber assortments [14][15][16]. Additionally, forest management considering the social cost of carbon favors increased tree size [17]. Since dominant individuals on productive rainforest sites can average diameter growth rates above 1 cm year −1 [18], trees can exceed small yarder payload capacities by age 40 and reach eight-wheel harvester limits by age 65 (Section 2.3). In northwestern United States rainforests, coast Douglas-fir (Pseudotsuga menziesii var menziesii (Mirb.) Franco) and western hemlock (Tsuga heterophylla (Raf.) Sarg.) trees are commonly harvested from coastal and inland mountain ranges at ages 40-50. Harvest at ages of 60-100 years is also of regular interest [19][20][21]. This combination of rapid growth, large trees, and rugged terrain creates need for steep slope harvest system models that remain accurate for tree diameters up to 100 cm or larger.
Because existing models of heavy harvest equipment apply primarily to trees of 60 cm diameter or smaller on slopes less than 50% (Section 2.1), this study models mechanized harvesting of larger trees on steeper slopes. We therefore (1) review and extend current equipment modeling approaches, (2) evaluate the proposed extensions by estimating uncertainty in predictions of harvest productivity and cost in two Douglas-fir and western hemlock case studies, (3) estimate harvest cost sensitivity to model input parameters, and (4) identify areas where additional research is needed to reduce uncertainty in Figure 1. Stump to mill flow of wood in three steep slope harvest systems. In system (a), a tethered harvester-forwarder pair provides cut-to-length harvesting suitable for thinning and continuous cover forestry, variable retention harvests, or clearfelling. Systems based on yarders (b,c) can be used for thinning but are most efficient in clearfelling [5].
Each piece of heavy equipment in a harvest system has size or weight restrictions on the wood it can handle (Section 2.1). Trees that are too large will not fit in felling or processing heads for cutting or may weigh too much for equipment booms to lift or drag. Logs too heavy for a forwarder's crane cannot be lifted into its bunk. Similarly, skyline deflection determines maximum payload weights in yarding. However, growing large trees on slopes of 50-100% is a frequent silvicultural objective within western North America's temperate rainforests for reasons such as conserving biodiversity, restoring wildlife habitat, carbon sequestration, and raising the financial value of timber assortments [14][15][16]. Additionally, forest management considering the social cost of carbon favors increased tree size [17]. Since dominant individuals on productive rainforest sites can average diameter growth rates above 1 cm year −1 [18], trees can exceed small yarder payload capacities by age 40 and reach eight-wheel harvester limits by age 65 (Section 2.3). In northwestern United States rainforests, coast Douglas-fir (Pseudotsuga menziesii var menziesii (Mirb.) Franco) and western hemlock (Tsuga heterophylla (Raf.) Sarg.) trees are commonly harvested from coastal and inland mountain ranges at ages 40-50. Harvest at ages of 60-100 years is also of regular interest [19][20][21]. This combination of rapid growth, large trees, and rugged terrain creates need for steep slope harvest system models that remain accurate for tree diameters up to 100 cm or larger.
Because existing models of heavy harvest equipment apply primarily to trees of 60 cm diameter or smaller on slopes less than 50% (Section 2.1), this study models mechanized harvesting of larger trees on steeper slopes. We therefore (1) review and extend current equipment modeling approaches, (2) evaluate the proposed extensions by estimating uncertainty in predictions of harvest productivity and cost in two Douglas-fir and western hemlock case studies, (3) estimate harvest cost sensitivity to model input parameters, and (4) identify areas where additional research is needed to reduce uncertainty in harvest system models. The following sections describe the cost and productivity modeling approach used for the three steep slope harvest systems in Figure 1, present harvest cost probability distributions for even-age coast Douglas-fir and western hemlock, and discuss limitations arising from the lack of large tree data and productivity models for some types of equipment.

Harvest Equipment Productivity Models
A literature search identified 20 sources (Table 1) with 44 productivity models of interest published in the last 25 years. Twenty-two of these models appeared to be errorfree and we were able to identify and correct errors in the equations of six others, resulting in Table 1. Eight of the remaining 16 models had errors that we were unable to repair and eight had forms unsuited to prediction. Therefore, we excluded these 16 models. We also excluded models lacking goodness of fit information, models with adjusted R 2 < 0.25, and sources older than 25 years. Among the resulting 20 sources, only the chainsaw productivity models included felling or bucking on slopes greater than 47%. No mechanized study indicated an average tree volume greater than 2.6 m 3 . Table 1. Cycle time definitions for productivity of harvest equipment used in steep slope mechanized logging and maximum slopes considered by previous studies. We propose two unified forms for cycle time models, one for machines with stemwise cycles and one for machines with roundtrip cycles (Equations (1) and (2)). loader stemwise no (1) pile one bunch of logs from processor (2) load one grapple onto truck none identified [45,49] log truck roundtrip no drive to unit, load logs into one or more bunks, drive to mill, unload logs [58] unclear but >11 [59,60] Log truck types considered in this study are mule trains with two bunks of logs, often ≤7.6 m and ≤10.1 m long, and long log trucks with one bunk of logs ≤ 12.5 m long. We also searched for sources on winch-assisted skidding but did not find published evidence of its uphill use on slopes above 35% [61,62].
Lack of productivity models covering slopes of 47-100% and trees with merchantable volumes up to 15 m 3 required these existing models be extended beyond their fitting ranges. We chose to extend cycle time models since most models that regressed productivity directly exhibited poor out of range behaviors, including prediction of negative m 3 PMh 0 −1 (merchantable cubic meters of wood per delay-free productive machine hour) or implausibly high productivities for their machine type. For machines with stemwise cycles (Table 1), the most common model forms for predicting the total cutting time T c for a given stem s were subsets of a per-stem equipment move time plus a power series of either the tree's merchantable volume V s or a highly correlated measure of tree size such as diameter at breast height. The per-stem move time is determined by the distance D move,s the equipment must move to reach stem s and the speed v s at which the equipment moves as a function of the local slope S s , traction T s available due to surface conditions, the amount of brush B present, and possibly other factors. As only one source [57] found evidence for higher-order terms in volume, we assumed the quadratic model form where b {0,1,2,3},m ≥ 0 are the felling and bucking coefficients for a particular machine of interest m and S t,m is the slope at which machine productivity begins to decrease. Models with quadratic terms commonly assume b 3,m = 0 but we allow b 3,m > 0 as in Hildt et al. [47] and use up to two quadratic terms (n q ≤ 2) to obtain productivity curves with shapes similar to those found by Visser and Spinelli [30] for trees of volume up to 5.2 m 3 . The slope term b 4,m (S s -S t,m ) is a linearization used in industry and is consistent with published estimates of slope effects [33,62].
For machines with roundtrip cycle times, a majority of sources ( Table 1) used models that implicitly or explicitly performed piecewise integration of mean unloaded and loaded travel speeds v s over a series of segments {s} with individual lengths D s . With the inclusion of time t spent loading and unloading, a model for roundtrip cycle time T c is While loaded and unloaded travel is frequently simplified to an average roundtrip speed over a single segment, more detailed approaches recognize uphill travel may be constrained by engine power P e and payload P l [58,62] as well as slope (S s ) and traction (T s ) [29,48,63]. The use of winch assist may increase [62] or decrease [36] travel speeds, potentially requiring tethered segments be modelled separately from untethered segments where winch assist is not used. Load and unload times are likely to increase with the number of timber assortments N s in the load [29] and the number of pieces (trees or logs) N p [4,47] transported, though the latter is routinely mitigated by grappling multiple pieces simultaneously or by yarding and swinging whole trees.
In addition to being specific to machine capabilities [28,29] model coefficients are likely to be tree species specific, both due to crown architecture [25,29] and differences in stem density and bark losses [64,65], which alter the amount of merchantable wood volume in a given payload weight. We provide specifics of the model parameterizations we use for feller-bunchers, harvesters, chainsaws, forwarders, yarders, processors, loaders, and log trucks with Douglas-fir and western hemlock in the Supplementary Material (Section S1). Since these parameters are extrapolated, we treat them as uncertain and include their uncertainty in estimates of harvest cost (Section 2.3).

Equipment Operating Cost and Harvest Contracting Business Models
Most harvest cost studies (Section 2.1) follow Miyata's method [66] for finding equipment operating costs. Under Miyata's method, costs are most sensitive to equipment depreciation, equipment utilization, and operator wages and compensation. We gathered a dataset of asking prices for harvest equipment offered for sale (n = 261) and, where sufficient data was available, estimated depreciation rates from nonlinear regressions of price versus operating hours. Machine utilization was estimated from figures reported in the literature review and operators' total compensation was estimated from job postings for positions in Oregon and Washington's temperate rainforests (a region of the northwestern United States where tethered equipment is used). We then extended Miyata's method to the consider the taxes, profit, and risk of a harvest contracting business beyond the cost of its equipment and operators. Similar to the productivity model parameterizations, we provide details of the business cost estimates for the 11 types of equipment we consider in the Supplementary Material (Section S2).

Harvest System Modeling and Uncertainty Estimation
A harvest system's financial cost is the total cost of its machines and their operators. In addition to the cost of operating the system to produce merchantable wood we included the costs C t for moving equipment into and out of a harvest unit, reopening roads, and other harvest related tasks (Supplementary Material, Section S3). When using cycle time productivity models (Section 2.1), a machine's wood production cost is found by multiplying the total cycle time T m over the n c cycles required to move the volume of wood harvested by the machine's cost per productive machine hour C m (US$ PMh 0 −1 in this study). The harvest cost for a system with n m machines is therefore Harvest costs are typically lowest when machines can operate independently of each other (also referred to as decoupled operation) as doing so maximizes their utilization and minimizes operating cost. We assume independent operation for feller-bunchers, harvesters, chainsaws, forwarders, and log trucks. Utilization of yarders, processors, and loaders, however, is likely interdependent. When one of these three machines restricted the productivity of the other two, we adjusted the costs of the affected machines accordingly (Section S3). When trees exceeded machine limits (Section S4), we included the cost of the chainsaw work needed to make wood accessible for mechanized movement.
We defined 72 productivity, 838 cost model, and 36 equipment limit parameters across the harvest systems considered (Sections S1-S4). We estimated harvest cost probability distributions empirically using 10,000 Monte Carlo draws of 68 parameters (48 productivity, 14 cost, and six limit). Given the lack of previous studies, we treated parameter ranges as uninformative Bayesian priors with uniform probability distributions estimated by Uusitalo et al.'s [67] expert judgement method. Simplification of Miyata models to uncertain operating costs per scheduled machine hour (US$ SMh −1 ) caused 94% of the reduction from the 946 parameters defined to the 68 Monte Carlo parameters (Section S2). The rest of the reduction was due to overlap among and quantization of forwarder parameters (Section S1.2) and holding biometric or harvest unit properties, such as bark losses and available soil traction, constant to avoid confounding harvest system variability with biological, seasonal, or unit-specific sources of variability.
Finally, we estimated cost model sensitivity to the 68 Monte Carlo parameters using total Sobol' indices [68] at N = 1500 (N(k + 2) = 105,000, Jansen estimator [69]) in order to prioritize model parameters for future investigation. We chose to use total Sobol' indices as they provide a comprehensive, variance-based measure of model sensitivity by considering a parameter's direct effect and its interactions with all other parameters.

Growth and Yield Modeling, Taper Equation Interactions, and Thinning Intensities
We initialized the northwest Oregon variant of the nonspatial, individual tree growth model Organon [70,71] with measurements from 30-year-old, even-age Douglas-fir and western hemlock permanent sample plots [18] and projected stand growth over the next 55-70 years. For each five-year timestep, we estimated per hectare harvest cost probability distributions for unthinned and thinned rotations by calculating stemwise cycle times for each individual tree present in the growth model and roundtrip cycle times for the volume of wood harvested per hectare. Stem cutting diameters and piece weights were estimated from using species specific taper equations, densities, and rates of bark loss [64,[72][73][74]. Financially preferred rotations, thinning intensities, and timings were identified by classic Faustmann land expectation values [75] calculated from mean timber prices over the past decade [76] net of harvest costs at a 4% discount rate [77]. We considered even-age rotations with up to three proportional thins and harvest intensities from 20-50% for each thin with step sizes of 1-3%.

Results
Physiological differences between coast Douglas-fir and western hemlock resulted in distinct silvicultural approaches, harvest cost distributions, and patterns of harvest equipment selection (Figures 2-5). Douglas-fir is a shade intermediate species capable of more rapid growth than the shade tolerant western hemlock (50-year site indices of 39.6 and 29.1 m, respectively, for the site used in this study [18]) has mean wood prices 22% higher [76], and log densities 9% lower [64] than hemlock. Douglas-fir rotations therefore supported up to three thins beginning at age 35, rather than at most one thinning for hemlock (Figures 2 and 4a,e). Douglas-fir trees grew to larger sizes and produced larger logs than western hemlock even though Douglas-fir rotations were 15 years shorter than hemlock rotations (Figures 3 and 5i-k). Harvest equipment therefore encountered larger stemwise quadratic cycle time contributions (Equation (1)) in Douglas-fir. Douglas-firs were also more likely to exceed mechanized equipment limits (Sections 2.3 and S4), resulting in greater costs for manual chainsaw labor beginning from age 45 ( Figure 6). Combined, these two effects led to well-defined minima in Douglas-fir harvest cost around age 50-60 whereas minima were absent or shallow for western hemlock (Figures 2 and 4a,b,e,f). Despite these differences, overall ranges of harvest costs were similar between the two species, though thinning in hemlock tended to be somewhat more expensive (Table 2).  [77]. We considered even-age rotations with up to three proportional thins and harvest intensities from 20-50% for each thin with step sizes of 1-3%.

Results
Physiological differences between coast Douglas-fir and western hemlock resulted in distinct silvicultural approaches, harvest cost distributions, and patterns of harvest equipment selection (Figures 2-5). Douglas-fir is a shade intermediate species capable of more rapid growth than the shade tolerant western hemlock (50-year site indices of 39.6 and 29.1 m, respectively, for the site used in this study [18]) has mean wood prices 22% higher [76], and log densities 9% lower [64] than hemlock. Douglas-fir rotations therefore supported up to three thins beginning at age 35, rather than at most one thinning for hemlock (Figures 2 and 4a,e). Douglas-fir trees grew to larger sizes and produced larger logs than western hemlock even though Douglas-fir rotations were 15 years shorter than hemlock rotations (Figures 3 and 5i,j,k). Harvest equipment therefore encountered larger stemwise quadratic cycle time contributions (Equation (1)) in Douglas-fir. Douglas-firs were also more likely to exceed mechanized equipment limits (Sections 2.3 and S4), resulting in greater costs for manual chainsaw labor beginning from age 45 ( Figure 6). Combined, these two effects led to well-defined minima in Douglas-fir harvest cost around age 50-60 whereas minima were absent or shallow for western hemlock (Figures 2 and 4a,b,e,f). Despite these differences, overall ranges of harvest costs were similar between the two species, though thinning in hemlock tended to be somewhat more expensive (Table 2).     Table S1).  Table S1). Table 2. Means and 95% intervals (in parenthesis) observed across 10,000 Monte Carlo parameter draws for thinning and regeneration harvests in Douglas-fir and western hemlock. Time totals are PMh 0 sums across all machines in the harvest systems considered (Figure 1)   Both extraction distance and slope significantly impacted harvest costs (p < 0.001). Increasing slope increased stemwise cycle times (Equation (1)), lowering productivity (Figures 3 and 5a,b,e,f) and, therefore, increased harvest cost as more machine hours were required to complete a harvest. Increasing slope from 65 to 80% decreased forwarder payload by 56% (Figures 2k and 4k, Equation (S1)), requiring 1.8 times as many forwarder turns to extract a given volume of wood. The corresponding increase in harvester hours was 1.13 times. Decreased payload was thus the primary cause of increased thinning costs on steeper slopes. Engine power did not limit forwarding speeds because driving speed while loading is more constraining than engine power. Similarly, extracting wood from 500 m long corridors increased forwarding times compared to 200 m long corridors, lowering forwarding productivity (Figures 3 and 5a,e) and also increasing thinning costs. Additionally, 500 m corridors exceeded the 280-350 m cable length limit of the add on winches available for use with eight-wheel harvesters and forwarders (Table S6). In such cases, harvesters and forwarders were supported by separate anchor machines, increasing costs compared to the use of an add-on winch powered by the harvester or forwarder.  The long log systems used with regeneration harvests were less sensitive to slope and distance than the cut-to-length, eight-wheel harvester-forwarder system we assumed for thinning. Since yarders, processors, and loaders were assumed to operate from landings or roadsides, their productivity was unaffected by the harvest unit's prevailing slope (Figures 3 and 5b-d,f-h). Similar to forwarding, yarding speeds did not appear to be constrained by engine power, possibly due to line speeds of partially suspended turns being limited by concern for hangups. Most cable yarders we surveyed could access the full length of 500 m corridors without cable extensions (Table S6). Additionally, the tracked feller-bunchers and harvesters considered were tethered by anchor machines and hence did require additional equipment to operate on corridors longer than 280-350 m. As slope and extraction distance increased, the long log systems exhibited an increasing cost advantage compared to a harvester-forwarder, cut-to-length system.
Yarder productivity was influential to harvest equipment selection in long log systems. Despite its 45% higher operating cost on an hourly basis (Table S4), most Monte Carlo parameter draws selected a swing yarder over an excavator-based yarder. Higher yarder productivity on 200 m corridors was associated with increased selection of a feller-buncher, even though the use of a feller-buncher typically increased chainsaw costs (Figures 2 and 4c,d,g,h and 7). The more rapid yarding allowed greater processor and loader utilization, making it more likely operating a feller-buncher with manual chainsaw bucking and a processor would be more cost effective than operating a single harvester at lower productivity (Figures 3 and 5b,f). With 500 m corridors, this effect reversed and harvesters were generally preferred. Because we assumed tracked harvesters were equipped with higher capacity processing heads than those currently available on an eight-wheel harvester, tracked harvesters were more likely to be preferred at increasing stand ages due to lower chainsaw costs (Figures 2, 4c,d,g,h and 6a,b,e,f). indicate which approach to staffing manual chainsaw work is most likely to result in the lowest cost for thinning (c,g) and regeneration harvests (d,h). In most cases, felling trees using mechanized equipment prior to manual bucking by a member of the cable yarding crew is preferred (Sections 1, 2.1, S2, and S4). Staffing bars do not reach a probability of one as chainsaw work is not always required. For example, no hemlock trees exceeded mechanized limits in thinning and regeneration harvests in hemlock did not exceed the assumed capabilities of a tracked harvester. Both extraction distance and slope significantly impacted harvest costs (p < 0.001). Increasing slope increased stemwise cycle times (Equation (1)), lowering productivity (Figures 3 and 5a,b,e,f) and, therefore, increased harvest cost as more machine hours were required to complete a harvest. Increasing slope from 65 to 80% decreased forwarder payload by 56% (Figures 2k and 4k, Equation (S1)), requiring 1.8 times as many forwarder turns to extract a given volume of wood. The corresponding increase in harvester hours was 1.13 times. Decreased payload was thus the primary cause of increased thinning costs on steeper slopes. Engine power did not limit forwarding speeds because driving speed while loading is more constraining than engine power. Similarly, extracting wood from 500 m long corridors increased forwarding times compared to 200 m long corridors, low- Figure 6. Distribution of harvest costs resulting from manual chainsaw bucking and felling of Douglas-fir (a,b) and western hemlock stems (e,f) exceeding the weight or diameter limits of mechanized equipment as a function of tree species and harvest equipment used. (c,d,g,h) indicate which approach to staffing manual chainsaw work is most likely to result in the lowest cost for thinning (c,g) and regeneration harvests (d,h). In most cases, felling trees using mechanized equipment prior to manual bucking by a member of the cable yarding crew is preferred (Sections 1, 2.1, S2 and S4). Staffing bars do not reach a probability of one as chainsaw work is not always required. For example, no hemlock trees exceeded mechanized limits in thinning and regeneration harvests in hemlock did not exceed the assumed capabilities of a tracked harvester. Sensitivity analysis of the parameters for the harvest cost model, other than slope and extraction distance, found that log transport (P l,kg , v, t roundtrip ) and equipment operating costs were consistently among the most influential parameters (Figure 7). Other influential parameters were the eight-wheel harvester constant and linear cycle time coefficients (b 0 , b 1 ) and the onset of slope induced decline in productivity (S t ). While eight-wheel harvesters were expected to be influential to thinning because other felling and processing machines were not considered, eight-wheel harvesters retained their influence in regeneration harvests where feller-bunchers, tracked harvesters, and processors were considered. Chainsaw parameters had limited influence because chainsaw costs were low in thinning and equipment selection for regeneration harvests limited chainsaw costs (Figures 2, 4c,d,g,h and 6b,f).  (2)). $ SMh −1 = operating cost per scheduled machine hour (Table S4). b0, b1, and St are stemwise cycle time coefficients (Equation (1)). Pl,kg = forwarding or yarding payload in kg (Section S1.2), , = power of forwarder payload in merchantable m 3 (Pl,c, Table S2). troundtrip = roundtrip trucking time (Table S3). Sensitivity analysis of the parameters for the harvest cost model, other than slope and extraction distance, found that log transport (P̅ l,kg, v̅ , troundtrip) and equipment operating costs were consistently among the most influential parameters (Figure 7). Other influential parameters were the eight-wheel harvester constant and linear cycle time coefficients (b0, b1) and the onset of slope induced decline in productivity (St). While eight-wheel harvesters were expected to be influential to thinning because other felling and processing machines were not considered, eight-wheel harvesters retained their influence in regeneration harvests where feller-bunchers, tracked harvesters, and processors were considered. Chainsaw parameters had limited influence because chainsaw costs were low in thinning and equipment selection for regeneration harvests limited chainsaw costs (Figures 2 and  4c,d,g,h and Figure 6b,f).

Discussions
The harvest cost model developed here attempts to provide reasonable, default parameterizations across a range of steep slope harvest systems. However, it relies on extrapolation from a limited body of results obtained from smaller trees on lower angle slopes (Section 2.1) and tethered harvesting equipment remains relatively new technology with some manufacturers, at least, making multiple changes per year [9]. Because all types of harvest equipment considered were sometimes selected, it appears likely which equipment is most preferred will depend on operating cost and productivity details of the specific crew members and harvest equipment available. It is also likely that the sizes, shapes, Figure 7. The 10 most influential harvest cost model parameters identified by Sobol' sensitivity analysis for Douglas-fir (a,b) and western hemlock (c,d) as a function of stand age and harvest type. Harvest unit corridor length and slope influence harvest costs as well but were excluded from this analysis as their effects are displayed in Figures 2-5. v = yarding line speed, both loaded and unloaded, or forwarding speed (Equation (2)). $ SMh −1 = operating cost per scheduled machine hour (Table S4). b 0 , b 1 , and S t are stemwise cycle time coefficients (Equation (1)). P l,kg = forwarding or yarding payload in kg (Section S1.2), P power l,c = power of forwarder payload in merchantable m 3 (P l,c , Table S2). t roundtrip = roundtrip trucking time (Table S3).

Discussions
The harvest cost model developed here attempts to provide reasonable, default parameterizations across a range of steep slope harvest systems. However, it relies on extrapolation from a limited body of results obtained from smaller trees on lower angle slopes (Section 2.1) and tethered harvesting equipment remains relatively new technology with some manufacturers, at least, making multiple changes per year [9]. Because all types of harvest equipment considered were sometimes selected, it appears likely which equipment is most preferred will depend on operating cost and productivity details of the specific crew members and harvest equipment available. It is also likely that the sizes, shapes, and slopes of harvest units will affect which harvest systems are preferred. In particular, the proposed model relies heavily on just two studies of yarders, one of excavator-based yarding on slopes less than 50% in Malaysia [51] and one of an older swing yarder operating in young Douglas-fir on slopes less than 60% [5]. In addition, the diversity of results for harvesters [26,28,[31][32][33][34]36,38] suggests no single model is likely to accurately represent all machines. Notably, increasing selection of tracked harvesters at later stand ages results from assuming they are fitted with larger capacity processing heads than eight-wheel harvesters (Table S6). With heads of equal capacity, a wheeled harvester may be preferred.
Given observed model sensitivities (Figures 2, 4c,d,g,h and 7), further study comparing productivity and operating costs of excavator-based and swing yarders appears likely to be effective in improving understanding of model uncertainty. Because their selection interacts with yarder productivity, comparative study of wheeled harvesters, tracked harvesters, and processors also seems likely to better quantify-and possibly reduce-model uncertainty. Understanding harvest contractors' reasons for selecting excavator-based yarders over other options may be valuable since, for example, the potential for more cost effective operation with a swing yarder is irrelevant if its more expensive purchase cannot be financed. While only grapple yarding in regeneration harvests was considered here, lateral yarding with chokers allows cable yarders to be used in thinning. While choker yarding is higher cost (Table S4) and lower productivity [5] than grapple yarding, it may provide more cost effective thinning than harvester-forwarder systems when slopes are steep enough and corridors long enough.
More complete understanding of tethered forwarder payloads also appears important. While forwarder tractive force was excluded from this study's Monte Carlo parameter draws and sensitivity analysis because it is primarily a property of the harvest unit, surface conditions can affect forwarder travel speeds by a factor of two [63] and payload changes of ±20% or more appear plausible. Since cut-to-length harvest systems can be cost competitive with long log systems (Figures 2 and 4a,b), it appears likely that steep slope conditions exist where harvester-forwarder pairs are more cost effective than cable yarder-based systems for regeneration harvests. Some harvester-forwarder operators suggest this is the case [39] and preferentially schedule steeper slope operations for seasons where traction is greater [9].
The apparent lack of loader models prevented assessment of harvest cost sensitivity to loader productivity parameters and may affect sensitivity to loader costs. Since loader operating costs are among the most influential parameters for long log harvest costs (Figure 7b,d) and would also occur when cut-to-length systems include loaders (Figure 1a) this suggests log loaders are an understudied type of harvest equipment. It follows that time and motion studies to develop loader cycle time models would be valuable to the completeness of harvest system models, creating an ability to quantify and prioritize loader uncertainties.
We caution that the results presented here do not consider site, machine, or operator specific uncertainty due to the confidential nature of private harvest contractors' operating costs. While contractors we work with have indicated consistency with their typical cut-to-length and cable yarding costs, they have declined to share details of their unitlevel production targets needed for more complete model calibration. Productivities reported in the literature (Table 1) sometimes differ by a factor of two or more, even for apparently comparable equipment and sites, and utilization can differ by 20% or more among crews [51]. Where possible, it is therefore recommended to use specific, rather than generalized, parameter estimates and to consider the parameterization details discussed in the Supplementary Material.
As with harvest cost models, growth model parameters are uncertain. However, parameter uncertainty has not been expressed for the growth model used in this study and we did not consider taper equation uncertainty. Tree growth and mortality is also subject to future climate uncertainty as well as episodic disturbance. Broader analysis to place harvest cost uncertainties into context with biotic and other abiotic uncertainties therefore appears desirable.

Conclusions
To our knowledge, this study is the first to propose a stump to mill harvest cost model enabling comparison of harvest systems across a range of tree sizes from 0-15 merchantable m 3 and slopes of 0-100%. While the model exhibits plausible behavior and we conducted both uncertainty and sensitivity analyses, it is subject to present limits in the knowledge of how tethered harvest equipment functions across tree size and slope. Therefore, the choices of model form and parameter values are somewhat speculative. Consequently, further study is needed to better understand the structure and parameterization of steep slope harvest cost models, particularly in regards to comparative interactions between harvest system types, harvest equipment types, specific machines, and operator abilities.