Influence of timber harvesting costs on the layout of cuttings and economic return in forest planning based on dynamic treatment units

Aim of study: To analyze the influence of harvesting costs on the distribution and type of cuttings when forest management planning is based on the dynamic treatment units (DTUs) approach. Area of study: A Mediterranean pine forest in Central Spain. Materials and methods: Airborne laser scanning data were used in area-based approach to predict stand attributes and delineate segments that were used as calculation units. Predicted stand attributes and existing models for diameter distribution and individualtree growth were used to simulate alternative management schedules for each segment for a 60-year planning horizon divided into three 20-year periods. Three alternative forest planning problems were formulated. They aimed to maximize or minimize net income, or maximize timber production with a constant flow of harvested timber. Spatial goals were used in all cases to enhance the clustering of treatments. Main results: Maxizing timber production without considering harvesting costs can be costly, even close to the plan that minimized net incomes. Maximizing net incomes led to frequent use of final felling instead of thinnings, placing cuttings near forest roads and creating more compact DTUs than obtained in the plan that maximized timber production. Research highlights: Compared to previous studies on DTUs, this study integrated felling and forwarding costs, which depended on distance to road and stand attributes, in the process of creating DTUs by means of spatial optimization. Additional keywords: spatial optimization; decision-making; forest economics; forest planning. Abbreviations used: ALS (Airborne Laser Scanning); CA (Cellular Automata); CC (Cut-Cut); CHM (Canopy Height Model); CNC (Cut-NonCut); DTM (Digital Terrain Model); DTU (Dynamic Treatment Units); FF (Final Felling); FNF (Final felling-No Final felling); SA (Simulated Annealing). Authors ́ contributions: Adrián Pascual conducted the analyses together with Timo Pukkala and Petteri Packalén, both involved in supervising the first author in the design of the study, data analysis and reporting. Annukka Pesonen was responsible for the segmentation and its documentation. Sergio de Miguel participated in the writing of the manuscript and supervision. All authors have read and approved the final manuscript. Citation: Pascual, A.; Pukkala, T.; de-Miguel, S.; Pesonen, A.; Packalen, P. (2018). Influence of timber harvesting costs on the layout of cuttings and economic return in forest planning based on dynamic treatment units. Forest Systems, Volume 27, Issue 1, e001. https://doi.org/10.5424/fs/2018271-11897 Received: 14 Jun 2017. Accepted: 07 Feb 2018. Copyright © 2018 INIA. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 International (CC-by 4.0) License. Funding: University of Eastern Finland; Academy of Finland (projects FORBIO-14970 and ADAPT-14907); European Union's Horizon 2020 MultiFUNGtionality Marie Skłodowska-Curie (IF-EF No 655815 to SdM). Competing interests: The authors have declared that no competing interests exist. Correspondence should be addressed to Adrián Pascual: adrian.pascual.arranz@uef.fi; adrian.pascual.arranz@gmail.com Forest Systems 27 (1), e001, 10 pages (2018) eISSN: 2171-9845 https://doi.org/10.5424/fs/2018271-11897 Instituto Nacional de Investigación y Tecnología Agraria y Alimentaria O. A., M. P. (INIA)


Introduction
Designing a forest plan is a complex task requiring the combination of many elements to achieve the management goals of the forest landowner (McDill, 2014). Contemporary forest planning methods are able to address a variety of simultaneous goals considering, for instance, timber harvesting targets and the spatio-temporal allocation of cuttings based on economic information (Pukkala, 2002). In these methods, preferences are converted into objectives and constraints when formulating multi-criteria decisionmaking problems (Diaz-Balteiro & Romero, 2008).
The first step of the planning process is to acquire information concerning the characteristics of a given forest management unit. For this purpose, airborne laser scanning (ALS) data are increasingly used since they represent a source of spatially-continuous information on forest attributes (McRoberts, 2006). Once presentstate information is available, alternative candidate solutions are compared to find the best possible management considering both preferences and consequences (Kangas & Kangas, 2002). Ensuring long-term sustainability of timber production has been a key element in the development of methods for forest management planning. In recent decades, the integration of additional criteria such as spatial objectives (Öhman, 2001) has increased the complexity of forest planning problems.
The integration of high-resolution spatial information from ALS techniques has eased and improved the inclusion of spatial goals in forest management planning (Baskent & Keles, 2005;Gobakken & Naesset, 2008). Moreover, forest inventory units can be delineated by means of ALS-based segmentation methods (Mustonen et al., 2008). This enhances the possibilities to use spatial optimization methods, for instance for delineating compact dynamic treatment units (DTUs; Holmgren & Thuresson, 1997; or improving the connectivity between old-growth forest patches (Tóth & McDill, 2008).
Although spatial goals have been used to cluster treatments when maximizing timber production (Öhman, 2001), it may be more meaningful to pursue spatial objectives indirectly, for instance via the effect of the location and size of harvest blocks on economic objectives (Bettinger et al., 2003). Previous research on spatial forest management optimization has addressed the implications of spatial planning constraints (Pukkala et al., 2014) in the economic performance of different problem formulations (Öhman & Eriksson, 2002;Augustynczik et al., 2016). However, previous studies involving DTU-based forest planning have ignored the dependence of harvesting costs on forwarding distance. Previous research in Spain developed functions for calculating the harvesting costs as a function of distance to road, slope of the terrain and forest attributes (Solano et al., 2007), which can be integrated in spatial forest planning studies (González-Olabarria & Pukkala, 2011).
The purpose of our study was to examine the influence of timber felling and forwarding costs on the location and type of cuttings. A forest management plan that maximized the difference between roadside value of harvested trees and harvesting cost was developed and compared to yield-maximizing plan, the latter representing the business as usual forest management in public forests of the study region. A third plan, minimizing net incomes, was also compiled as the worst possible scenario to properly locate the yield-oriented plan within the range of variation in the economic performance of alternative forest management options. The plans were compared according to their economic performance, type and intensity of cuttings, and the clustering of treatments.

Study area and field data
The study area of 1,059 hectares is located in the mountainous forests of the Iberian System, in Castilla y León (Spain) (41°51´N, 3°15`E). It contains a mixture of dense and sparse pine stands where Pinus nigra Arn. is the main species. The slope of the terrain is mostly between 5 and 15%. About 1% of the area is steeper than 35%, which is regarded as the limit for timber harvesting.
Forest management is oriented towards sustained production of timber giving increasing importance to commercial thinning. This is partly due to the increasing value of thinner timber assortments. Usually, an early-stage thinning is prescribed 20 years after the establishment of the natural regeneration followed by one or two thinnings. Then, seed tree cut is prescribed. The rotation length for naturally regenerated P. nigra forests is between 90 and 130 years depending on site productivity.
Field data were collected using a systematic grid of 116 circular plots (12.6 m radius). The required forest attributes needed for the planning system were number of stems per hectare, stand basal area, stand age and dominant height. Stand age was estimated using tree ring analysis on the tallest tree within each plot, and the height of the same tree was taken as the dominant height of the plot (Table 1). Forest road information was provided by the Forest Service of Soria province. It included data on 17 well-conserved forest roads (18,665 m in total, 17.6 m/ha) surrounding and intersecting the study area.

ALS data and segmentation
At the same time as the field work was carried out, ALS data were collected using the Leica ALS60 II laser scanning system. The average pulse density was 2 pulses/m 2 . The scanner recorded up to 4 echoes per pulse. The digital terrain model (DTM) was interpolated at a 1 m 2 pixel size from echoes classified as ground echoes using the method proposed by Axelsson (2000). The elevation of the raster cell was the average elevation of those ground echoes that were located within the cell. The canopy height model (CHM), also with the pixel size of 1 m 2 , was interpolated by searching the highest ALS echo within a radius of 1.6 m (Euclidean distance) from the centre of each pixel. Then, the DTM was subtracted from the ALS echoes to convert absolute heights to heights above ground. Pixels that remained empty were filled by calculating the average of the highest echoes from non-empty neighbouring pixels. The resulting CHM was used for delineating homogeneous segments in terms of stand structure (Fig.1). Forest roads were considered as 7-m wide polygons in segmentation (see right panel in Fig. 1) and the bordermost forest segments started exactly from the side of roads. Overall, 2,947 segments (average size 0.36 ha) were used as forest inventory units. The reader is referred to Pascual et al. (2017) for a detailed explanation on how the segmentation was done and how ALS metrics were selected and forest attributes predicted for the segments.

Felling and forwarding cost functions
We used the cost functions (Eqs. [1]-[2]) presented in earlier studies under similar conditions (Solano et al., 2007;González-Olabarria & Pukkala, 2011) to calculate the felling costs (including debranching and crosscutting) and forwarding costs to the nearest forest road. Slope, distance to road and the diameter of the harvested tree were the predictors of the cost functions (Fig. 2).

C felling
[1] where C felling is felling cost (€/m 3 ), C forwarding is forwarding cost (€/m 3 ), d is dbh of the harvested tree (cm), Slope is the average slope of the stand (%), and Distance is distance to road (m). The distance from a segment to the nearest forest road was computed along the surface of the terrain considering the slope. The slope was calculated from ALS data as the mean slope of the segment.

Planning problems and spatial optimization
Since the predicted forest attributes consisted of segment-level information, we used diameter distribution models to obtain individual-tree attributes. Then, the individual-tree models presented in Pascual et al. (2017) were used to predict stand dynamics along a 60-year planning horizon divided into three 20-year periods. For the 2,947 calculation units, we simulated alternative nearly optimal even-aged management schedules using the Monte forest vers. 6.0 planning software (Palahi et al., 2004).
To create the instructions for nearly optimal management, stand dynamics were used to predict the value increment of a large set of segments. Growing stock volume and value, stand basal area, site index and stand age were also calculated for each segment. Then, a regression model showing the value increment as a function of the above-listed stand variables was fitted and used to develop instructions for economically optimal management. The instructions showed the basal area at which thinning should be conducted and the mean diameter at which final felling should be done. The stand was considered financially mature for cutting when its relative value increment fell below 2%.
The simulation used these instructions for rotation length and thinning basal area. Seed tree cut followed by the removal of seed trees during the next period was simulated when the mean diameter was larger than the instruction while thinnings were simulated when stand basal area was higher than the thinning basal area threshold of the instruction. Then, rotation diameter and thinning basal area were multiplied with constants slightly larger or smaller than 1, and the simulation was repeated to obtain several treatment schedules per segment. Always when basal area exceeded the thinning limit, three thinning intensities were simulated: light (20% reduction in stand basal area terms), normal (30%) and heavy thinning (40%). This also increased the number of alternative treatment schedules. The 2,947 segments of the case study area had 60,581 different simulated schedules, which is 21 schedules per segment on average.
Three alternative forest planning problems were formulated and solved to examine the effect of harvesting costs on forest plan: i) maximize timber production (MaxPro); ii) maximize net income (MaxNetInc) and iii) minimize net income (MinNetInc). Since MaxPro ignored harvesting costs, cuttings were allocated without considering the effect of slope and forwarding distance on harvesting costs. MaxNetInc, which aimed at maximizing the total net income (roadside price of harvested trees -harvesting costs), took these factors into account. The MinNetInc formulation was not realistic regarding practical forestry, but it showed the lower bound of net income and making it possible to evaluate the MaxPro plan from the net income point of view.
Eight management objectives were included in each forest planning formulation. The principal management objective is referred to with symbol M as follows: • MaxPro: Maximization of standing volume at the end of the plan (M 1 ).
• MaxNetInc: Maximization of net income during the plan (M 2 ).
• MinNetInc: Minimization of net income during the plan (M 3 ).
The three formulations shared the following seven non-spatial and spatial management goals: • Removal of 50,000 m³ during each 20-year period (R 1 , R 2 , R 3 ).
• Maximization of the proportion of cut-cut borders of adjacent segments (CC).
• Maximization of the proportion of cut-cut borders of adjacent segments prescribed as final fellings (FF). • • Constant timber flow along the plan was required (objectives R 1 , R 2 and R 3 ) to ensure sustained timber production and steady economic return. Preliminarily, we tested three possible harvesting scenarios (25,000 m 3 , 50,000 m 3 and 100,000 m 3 ). The results revealed that the difference between initial standing volume (253, 016 m 3 ) was very close to the ending volume when 50,000 m 3 was harvested during every 20-year period (ending volume was 256,043 m 3 ). Therefore, we set the harvesting target as 50,000 m 3 in all problems.
The aggregation of harvest blocks was pursued via four simultaneous spatial goals that promoted the clustering of segments having a similar cutting prescription. CC and CNC aimed at aggregating all cuttings while FF and FNF specifically aggregated final fellings to avoid isolated segments prescribed as final felling within a thinning block. The aggregation was based on maximizing the shared border of adjacent cut-cut segments (CC and FF). In addition, the proportion of cut-uncut boundary was minimized to create compact and round-shaped treatment units (objectives CNC and FNF).
The cellular automaton (CA) of  was used as the optimization method since previous research has shown the good performance of decentralized heuristic methods in large spatial problems (Pukkala et al., 2009) and, specifically, to create compact harvest blocks (Pascual et al., 2017). This two-phase CA included segment-level and neighbourhood-related goals (local function) and forest management unit level goals (global function). The local function included the spatial Minimization of the proportion of cut-non-cut borders of adjacent segments prescribed as final fellings (FNF).
Minimization of the proportion of cut-non-cut borders of adjacent segments (CNC).
objectives (CC, FF, CNC and FNF) and the primary objective (M) of the respective formulation. The CA calculates the probability of mutation and innovation as a function of iteration number. Firstly, one of the simulated treatment schedules is randomly selected for each segment. Then, a random number from 0 to 1 is drawn from a uniform distribution. If the random number is smaller than the current mutation probability, a mutation takes place . Innovations are assigned in the same way. If innovation occurs, the management schedule of the segment is replaced by the one which maximizes the following local function (Eq. [3]): where U is the objective function value, s (sign) is +1 for MaxPro and MaxNetInc plans and -1 for MinNetInc plan and Mmax is the largest per hectare value of the principal objective variable among all schedules of all segments. Preliminary tests were used to assign such weights for the spatial objetives that ensured a good harvest block layout. In this regard, the weights of CNC and FNF were progressively increased until compact harvest blocks were obtained.
In the second phase of the CA, a global objective function (Eq. [4]) was added to the local one, expressing the cutting targets (R 1 , R 2 and R 3 ). The principal objective variable was included also in the global objective function as this resulted in slightly better solutions than having the non-spatial objective only in the local objective function: where p 1 , p 2 , p 3 and p 4 are the sub-priority functions for the cutting targets and M (Fig. 3).
In the global optimization phase, treatment schedules were re-assigned until the cuttings targets were met. The following objective function, including both the local and the global function, was used in the second phase of the CA run: where OF is the objective function value, a is the area of the segment, A is the total area of all segments, and b is the weight of the global priority function P. Once the local optimization phase was completed, the global phase started with Eq.
[5] as the objective function, so that the value of b was increased from zero by 0.01 in every iteration until the global utility (Eq. [4]) reached a value, which could no longer be improved. At this point, the even-flow harvesting targets were met.

Assessment of forest management plans
The accumulated net income during the planning period was used as an indicator to compare the three plans in financial terms while timber production was assessed based on the standing growing stock volume at the end of the 60-year plan assuring that the 50,000 m 3 harvesting targets were fulfilled for each 20-year period. The spatial distribution of prescriptions was analysed considering whether cuttings were final fellings and overstory removals, or thinnings. The prescribed area for each treatment and the spatial layout of the resulting DTUs were compared to identify changes in the location of prescriptions and their proximity to forest roads. The mean diameter of harvested trees in the resulting prescriptions was assessed together with the proximity of prescriptions to forest roads.

Management objectives
As expected, the MaxNetInc plan resulted in the highest total net income (6.24 million €). For MaxPro and MinNetInc plans, the net incomes were almost the same: 5.67 and 5.56 million €, respectively. In terms The analysis revealed how the MaxNetInc plan progressively diverged from the other plans from the first to the third period (Fig. 4).
The total amount of treated area needed to fulfil the cutting requirements in income-oriented planning (MaxNetInc) was 964.2 ha, 7% more than in MaxPro (898.2 ha) and 30% less than MinNetInc (1,253.2 ha). Differences in cutting type were largest during the second and third periods ( Table 2). During the first 20year period, the MaxNetInc and MaxPro plans prescribed many heavy thinnings whereas the MinNetInc plan employed light thinning over large areas. The cuttings of MaxNetInc were more concentrated in dense and mature forests during all periods compared to MaxPro, but especially during the second and third period. Fewer segments were treated with final felling in MaxPro, leading to increased thinning area. The mean diameter of harvested trees increased toward latter 20year periods in all cases, reaching up to 45.7 cm in the MaxNetInc plan, 6.3 cm more than in MaxPro plan.

Spatial layout of DTUs
The location of prescriptions was visually assessed by displaying the resulting DTUs of the three plans (Fig.  5). Compared to the MaxPro plan, prescriptions of the MaxNetInc plan tended to be more clustered along the road network during the second and third periods. This pattern can be seen for instance in the south-eastern part of the study area, where the location and type of prescriptions varied considerably among the plans. In the second period, more segments were prescribed as seed tree cut in MaxNetInc than in MaxPro plan, as incomeoriented plan preferred cuttings where the harvesting cost per cubic meter was low (harvested trees were large). The spatial layout of the harvest blocks showed an increasing level of compactness in MaxNetInc, while, in MaxPro, prescriptions were more scattered and less compact than in MaxNetInc and MinNetInc.

Forwarding distance and size of harvested trees
The observed differences in the allocation of prescriptions were further assessed by comparing the mean diameter of harvested trees and forwarding distance (Fig. 6). The mean distance to road of all cuttings was 253.9 m in MaxPro, which is clearly more than in MaxNetInc (168.8 m). In the first period, the scatter plots of harvested tree size and forwarding distance were very similar in MaxNetInc and MaxPro plans (Figs. 6a,b). The size of the harvested trees considerably increased in MaxNetInc during the second and third period (Figs. 6c,d,e,f). In the second period, cuttings prescribed in MaxNetInc were located close to forest roads (600 m maximum), while in MaxPro cuttings were prescribed up to 1,300 m from the nearest road.

Discussion
Three patterns emerged when net incomes were maximized (MaxNetInc) instead of timber production (MaxPro). In MaxNetInc plan, i) final fellings were prescribed more often, ii) the average size of harvested trees was larger and iii) cuttings were located more in the proximity of roads.
The observed differences between plans were more pronounced during the second and third 20-year period, which is in line with previous research (Öhman & Lämås, 2003). The lack of large differences in the first period between MaxNetInc and MaxPro can be partly explained by the design of the existing forest road network, which has been built to provide access to currently mature forests. After the best stands near roads had been harvested in the first period, the differences between MaxPro and MaxNetInc plans started to increase.
The spatial layout of the presented DTUs confirmed the previous results as prescriptions tended to follow forest roads, boosting the compactness of the harvest blocks of the MaxNetInc plan. The presence of isolated harvest blocks was lower in MaxNetInc than in MaxPro. Although the clustering of harvests was pursued in all three management plans, the integration of harvesting cost equations contributed to a better clustering of

MaxPro
First Period treatment units. Taking into account that maximization of net incomes tended to increase the area of final fellings, additional constraints may be considered to prevent the DTUs from becoming too large and not conforming with country-specific forest management regulations (Tóth & McDill, 2008;Carvajal et al., 2013). Large continuous areas of final fellings may also be detrimental for other ecosystem services, for instance recreational values and maintenance of biological diversity. In our study, we used a decentralized optimization algorithm to overcome the computational complexity of the planning problem in which the number of calculation units and treatment schedules was high (Pukkala et al., Figure 6. Distribution of mean diameter of harvested trees and distance to road for MaxNetInc during the first (a), second (c) and third period (e); and for MaxPro (b, d, f). The green (+) and red (x) marks show the mean distance to road and the overall mean diameter of trees in the prescribed cuttings for MaxPro (+) and MaxNetInc (x). 2009). High number of calculation units often prevents the use of exact methods such as integer programming. Before implementing the CA-based optimization, we tested simulated annealing algorithm (SA) as previous research has reported its good performance in similar problems (Pascual et al., 2016;Borges et al., 2014;Jin et al., 2016). The SA results were rather similar to those obtained with CA only for the MaxPro planning problem. In the two other problems (MaxNetInc and MinNetInc), harvest blocks were not as well clustered as with CA and, consequently, it was more difficult to see the effect of harvesting costs on the forest planning layout of cuttings when using SA. Our results confirmed those from previous research in relation to the better performance of CA compared to SA in spatially-explicit forest planning when the number of calculation units is high Mathey et al., 2007).
Our ALS dataset played a major role in our study as it was used in the delineation of forest inventory units by means of segmentation and for modelling and prediction of the initial values of forest attributes. The economic feasibility of ALS inventory depends on many factors, such as forest structure and information needs. The target area must also be large enough because a set of local field sample plots is needed and the required number of sample plots is more or less fixed although the inventory is enlarged substantially. If stand attributes are needed by tree species the required number of sample plots is many fold. ALS data acquisition is also expensive in small areas, but sometimes ALS data collected in other campaigns, e.g. for nationwide DTM purposes, can be used. ALS inventory is already considered to be economically viable both in plantations and semi-natural forest (see e.g. chapters 11-13 in Maltamo et al., 2014). The delineation of forest inventory units by means of segmentation can be done without any field data but without predicted forest attributes its value is limited.
By using full-coverage ALS data, we were able to estimate the spatial variation in harvesting costs across the study area. Previous research on spatial forest planning problems have relied on indirect methods such as external data to assign harvesting costs to each calculation unit (Öhman et al., 2011), or by assuming that machinery displacements are closely related to harvesting cost variation (Augustynczik et al., 2016). Variation in the quality of forest road network was not an issue in our study. This variation may affect the cost of transporting timber from roadside to saw and paper mills, which in turn would affect cuttings so that more wood should be forwarded to good-quality roads.
The economic return of the MaxPro plan was very close to that of MinNetInc, the worst possible scenario in economic terms. This means that forest management for maximal timber production was very costly.
Discounting, which was not done in our study, would affect the results so that cuttings are scheduled earlier when discount rate increases.
The differences between MaxPro and MaxNetInc plans showed that the degree of clustering between forest inventory units composing a harvest block may partly depend on management objectives which do not directly aim at aggregating treatments. Even when the purpose of the management is not entirely incomeoriented, the delineation of treatment areas should be flexible, responding to stand dynamics and consequent changes in stand structure, and combining multiple objectives. In this regard, the effect of maximising the combined production of timber, biomass and non-wood products (e.g., edible mushroom) on treatment decisions should be addressed in further studies. Consideration of harvesting costs in problem formulation (Solano et al., 2007) contribute to not only management objectives, but it also reinforces the clustering of management prescriptions, reducing the fragmentation of the forest.