Two 40 m by 40 m sample plots, one in the Eno village of Joensuu and the other in Revonkylä (also Joensuu) were used in this study (Table 1, Fig. 1). The plots belong to the research plots of the KESTO project, which investigates the implementation alternatives and harvesting costs of CCF, as well as soil compaction and logging damages to remaining trees. The KESTO plots are 20 m wide. When the plots are harvested, a strip road is opened in the middle of the plot. The two 40 m 40 m plots used in this study were obtained by joining two adjacent 20-m-wide plots. The EUREF-FIN-TM35FIN coordinates of the plots were: Eno 6964513 (N), 661689 (E); Revonkylä 6953432 (N), 678946 (E).
The plots were measured in stands where the landowner has decided to use CCF. The plots have already been thinned but this study used pre-thinning data. All trees taller than 1.4 m were measured before the thinning treatment for species, dbh, and location (Table 1, Fig. 1). Tree height was measured for every fourth tree. The logging injuries of all remaining trees were assessed after the thinning treatment. The Cartesian coordinates of the trees were determined by a distance from a reference point and an angle from a reference direction (north). For this study, the distances and angles were converted into x and y coordinates.
Table 1
Species composition of the sample plots of Eno and Revonkylä. N = number of trees per hectare, G = basal area (m2ha− 1), D = basal-area-weighted mean diameter (cm).
|
Eno
|
Revonkylä
|
|
N
|
G
|
D
|
N
|
G
|
D
|
Scots pine
|
6
|
0.3
|
24.2
|
150
|
4.7
|
22.5
|
Norway spruce
|
1294
|
22.4
|
24.9
|
2888
|
10.4
|
13.9
|
Silver birch
|
31
|
1.1
|
30.1
|
0
|
0
|
-
|
Downy birch
|
269
|
5.4
|
20.2
|
806
|
8.8
|
17.1
|
Aspen
|
100
|
5.9
|
31.8
|
0
|
0
|
-
|
Alder
|
25
|
0.2
|
15.3
|
0
|
0
|
-
|
Rowan
|
0
|
0
|
-
|
19
|
< 0.1
|
1.9
|
Total
|
1725
|
36.2
|
25.6
|
3869
|
24.0
|
16.7
|
Table 1, Fig. 1
Methods
The method developed and tested in this study was based on the alternated use of two optimization methods (Fig. 2). The higher-level method optimized the cutting years and the thinning intensity curves of trees for which the harvest decision was not optimized individually. The formula of the thinning intensity curve was (Sun et al. 2023):
$$TI\left(d\right)= \frac{1}{{\left[1+{a}_{1}\text{e}\text{x}\text{p}\left(-{a}_{2}(d-{a}_{3})\right)\right]}^{\left(\frac{1}{{a}_{1}}\right)}}$$
1
where TI(d) is the proportion of removed trees for breast height diameter (dbh) d, and a1, a2 and a3 are parameters. For each thinning treatment, the higher-level method optimized the following four variables: number of years from the beginning or previous cutting, and parameters a1, a2 and a3 of the thinning intensity curve. Therefore, the number of variables optimized at the higher level was equal to four times the number of optimized cuttings. The algorithm used at the higher level was differential evolution (Storn and Price 1997; Pukkala 2009; Jin et al. 2018).
Figure 2
Combinatorial optimization was used for trees for which the cutting event was optimized individually. For example, if the number of optimized cuttings was 3, each tree had four possible decisions: (1) removed in the first cutting, (2) removed in the second cutting, (3) removed in the third cutting and (4) not removed in any of these three cuttings. Finding the optimal combination of these four options for all trees is a regular combinatorial problem, for which many methods are available. We used simulated annealing (e.g., Bettinger et al. 2002), but integer programming (Pascual 2021), genetic algorithm (Fransson et al. 2020), tabu search, and cellular automata (Packalen et al. 2020), among others, could be used as well.
Two criteria were used to select the trees for tree-level optimization: initial breast height diameter (dbh) of the tree, and number of cuttings in which individual-tree optimization was used. For example, it was possible to optimize the cutting decision at the tree level only for trees larger than 20 cm in dbh, and only for the first cutting. In this case, Optimized thinning intensity curves would be used for smaller trees in the first cutting and for all trees in later cuttings.
Combinatorial optimization formed the lower level of the hybrid method (Fig. 2). Each combination of cutting years and thinning intensity curves that were tested by the higher-level algorithm was passed to the lower level. The lower-level algorithm allocated the trees selected for tree-level optimization to different cutting events (including the no-cut option).
Each combination inspected at the lower level was simulated, taking the cutting years and the thinning rates of the small trees from the higher-level solution. This simulation computed the objective function value of the tested management schedule. When net present value was maximized, the net income of a cutting event was calculated as the difference between the roadside value of the harvested trees and harvesting costs. Harvesting costs were calculated by using time consumption functions for the harvester and forwarder (Rummukainen et al. 1995).
The functions of Rummukainen et al. (1995) we used to calculate the harvester time of each harvested tree as a function of tree size (stem volume) and cutting type (clear felling or partial cutting). The time consumption of the forwarder depends on the total removal (m3/ha, converted into cubic meters per 100 meters of strip road, assuming a 20-m distance between strip roads), cutting type (partial cutting or clear-felling), terrain class (assumed normal) and forwarding distance (200 m assumed). The harvester and forwarded times were multiplied by the hourly cost of the machine (100 €/h for the harvester and 70 €/h for the forwarder). The road size values were 65 €/m3 for conifer saw logs, 50 €/m3 for birch saw logs, and 30 € for conifer and birch pulpwood. Aspens were crosscut to pulpwood logs, the roadside of which was 25 €/m3. The third assortment was energy wood with a roadside price of 20 €/m3 for conifer and birch, and 15 €/m3 for aspen.
The taper models of Laasasenaho (1982) were used to calculate assortment volumes. The over-bark minimum top diameter of the saw log was 15 cm for pine, 16 cm for spruce, and 18 cm for birch. The top diameter was 8 cm for pulpwood logs and 3 cm for energy wood logs. The minimum piece length was 4.3 m for conifer saw logs, 3.4 m for birch saw logs, and 3 m for pulpwood and energy wood logs.
The smallest accepted harvest removal was 50 m3ha− 1 per thinning since timber buyers may not be interested in buying the timber if the removal is smaller. The maximum removal was 150 m3ha− 1 because very high removals make the remaining stand vulnerable for instance to wind throws and may cause much damage to advance regeneration. The lowest accepted post-cutting basal area was 8 m2ha− 1, to guarantee that the optimized cutting schedule represented CCF.
The NPV of the ending growing stock (remaining growing stock after simulating the last cutting) was calculated with the formula of Pukkala (2022):
NPV end = exp(b1 + b2ln(TS) + b3Peat + b4ln(r)ln(TS) + b5ln(D×G×LogPrice)) (2)
where NPV end is net present value (€ ha− 1), TS is temperature sum (degree days), Peat is indicator variable for peatland forest (1 for peatland, 0 otherwise), r is discount rate (%), D is basal-area-weighted mean diameter of trees (cm), G is stand basal area (m2ha− 1), and LogPrice is the roadside price of saw log (€ m− 3). The parameters (b1–b5) are different for different size fertility categories (herb-rich, mesic, sub-xeric, xeric, barren).
The growth models of Pukkala et al. (2021) were used to simulate the diameter increment and survival of the trees, and the ingrowth of the stand. The time step of the simulation was 5 years, corresponding to the time step of the models. The simulation was conducted at the tree level, i.e., diameter increment and survival of each tree of the plot were simulated individually. The ingrowth model predicted the number of trees per hectare that will pass the 1.3 height during the next five years. This prediction was converted to the number of trees in a 40 m by 40 m plot (0.16 ha). For example, 50 ingrowth trees per hectare correspond to 8 new trees in a 0.16-ha plot. These new trees were placed in the plot in random locations (the x and y coordinates were drawn from uniform distribution U[0,40]).
The survival model is a logistic model showing the probability that a tree survives for 5 years. Survival was simulated by comparing the survival probability to a uniform random number (U[0,1]). Simulating survival in this way is stochastic, which is not ideal in optimization because the same set of optimized variables may result in different NPVs in different simulations. This would make it more difficult for the optimization algorithm to locate the optimum. To avoid this problem, a uniformly distributed random number (U[0,1]) was generated for each tree at the beginning of the optimization and this number was kept unchanged for the whole search process. The survival probability was compared to this random number.
The thinning intensity (Eq. 1), when applied at the individual tree level, is also a probability (probability of removal). Because of this, another uniformly distributed random number was generated for each tree at the beginning of the simulation optimization process. This number was compared to the value obtained from the thinning intensity curve to decide whether the trees should be removed.
All ingrowth trees also obtained the two tree-specific uniform random numbers. To guarantee deterministic simulation results for ingrowth, a set of random numbers for possible ingrowth trees was produced at the beginning of the optimization. Therefore, a certain ingrowth tree always obtained the same random numbers in repeated simulations of the same cutting schedule.
In most analyses of this study, NPV was maximized with a 3% discount rate. However, since one motivation for using CCF is the intention to increase the diversity and resilience of the forest, multi-objective optimizations were also conducted where the management objectives were NPV, tree species diversity, and the diversity of tree sizes. The following utility function was maximized:
U = w1(NPV/NPVmax) + w3(Shannon/Shannonmax) + w3(ModGini/ModGinimax) (3)
where NPV, Shannon, and ModGini are, respectively, the NPV, Shannon index and modified Gini index of the management schedule and NPVmax, Shannonmax, and ModGinimax are the maximum possible values of these variables in the stand. The maximum values were obtained from single-objective optimizations. Shannon and ModGini indices were the average values of these indices during the optimized period. The modified Gini index was obtained by multiplying the Gini index calculated from tree diameters (Valbuena et al. 2012) by the range of diameters (maximum dbh – minimum dbh). This was because the Gini index only describes the evenness of tree size and ignores the range of dbh variation (Valbuena et al. 2012).
The Shannon index (S) was computed from the basal areas of five different species or species groups:
$$S= -\sum _{i=1}^{5}{p}_{i}\text{l}\text{n}\left({p}_{i}\right)$$
4
where pi is the proportion of species i of the stand basal area. The subscript i refers to different species as follows: 1 pine, 2 spruce, 3 silver birch, 4 downy birch, and 5 other species (aspen, alder, or rowan).