Optimizing Thinnings and Rotation of Scots Pine and Norway Spruce Mixtures

The study describes a simulation-optimization system which uses spatial models for diameter and height growth, crown ratio and tree mortality for Scots pine and Norway spruce mixtures. The optimal oneand two-thinning regimes of six initial stands with varing species composition were solved by using nonlinear optimization. The soil expectation value (SEV) at 3 % interest rate was used as a management objective. The regimes are determined by taking into account the stand basal areas before the thinnings, the removal percentages for small, medium-sized and large pines and spruces, and the stand basal area before the final felling. The greatest SEV (8900 FIM ha –1) was attained with the initial stand where the proportion of pines was 65 % of the number of the stems. In the two-thinning regime, the first thinning was conducted at the age of 39 years when the stand basal area was 37 m 2ha–1 and the dominant height was about 15 m. After the thinning, the basal area was 27 m 2ha–1. Spruces were thinned from below, but both small and large pines were removed. The second thinning was 8 years later and much heavier: the stand basal area was decreased from 35 m 2ha–1 to 18 m2ha–1 by removing both small and large pines and spruces. When the optimal two-thinning regime was compared to the regime presented by Forest Centre Tapio, the loss of SEV was about 30 % (6070 FIM ha–1) in the case of thinnings from below, and about 20 % (7250 FIM ha –1) in the case of thinnings from above.


Introduction
In Finland, 65 % of the total land area are classified as forest land (Metsätilastollinen vuosikirja 1997), amounting in about 4 hectares of forest per each Finnish citizen.Forestry is a very important source of income for the Finnish national economy.One way of improving the growth and economic return of the forests on the sites of medium fertility is to grow mixed forests, e.g.mixed forests of Norway spruce (Picea abies (L.) Karst.) and Scots pine (Pinus sylvestris L.) (Jonsson 1962, Poleno 1981, Pukkala et al. 1994, 1998b, Vettenranta 1996) or mixed forests of birch (Betula sp.) and conifers (Mielikäinen 1980(Mielikäinen ,1985)).
There are several reasons why growing mixed forests can be more profitable to grow than growing pure stands (see e.g.Kelty et al. 1992).The species in mixed forests may have different requirements of light and nutrients and a different influence on the soil which result in the mixed forest effect (Mikola 1954, Sepponen et al. 1979and Kelty 1992).The different growth rhythm during a rotation also leads to the better productivity if the proportion of species is adjusted reasonable (Pukkala et al. 1994, Vettenranta 1996).Mixed forests also have a greater selection of treatment alternatives than pure stands, and the variation in prices of timber assortments brings along new possibilities to increase economic return.
The species composition has both short-and long-term effects on the growth dynamics as well as on the profitability of a certain treatment decision.The litter and the micro-climate affect the soil condition in the long run.Thus, the species adjustment may affect the regeneration and the growth rate of the next tree generation.Some treatments may have quite a short-term effects (e.g.intensity of thinning), while others (e.g.clearing and thinning of a sapling forest) may have crucial effects on the stand development even over rotation.
For reasons mentioned earlier, discovering the optimal treatment for a mixed stand is a more complicated task than that for a pure stand.The effects of all possible treatment decisions in mixed forests are almost impossible to detect with empirical studies, due to the numerous treat-ment alternatives and the great variation in the structure of mixed forests.Thus, the most profitable thinning regime and the structure of an initial stand can be determined only by simulation and optimization.
In the previous studies on the optimal thinning regime of coniferous mixtures, the thinning from above has been found to be more profitable thinning type than the thinning from below (Vettenranta 1996, Pukkala and Miina 1997and Pukkala et al. 1998b).However, the knowledge of the post-thinning development of a stand is uncertain, especially when a dense stand is thinned from above, due to the less vigorous vitality of the remaining earlier suppressed trees.In growth models, the thinning reaction of remaining trees have been predicted by e.g. a crown ratio (Mielikäinen 1985, Hynynen 1995, Vettenranta 1999), ad hoc-models (Salminen 1995), thinning response variable (Hynynen 1995, Jonsson 1995) and harvested competition (Mielikäinen 1978, Bucht 1981, Pukkala 1989, Tham 1989, Pukkala et al. 1998a).
In the growth models by Vettenranta (1999) in which the vitality of trees has been taken into account with the aid of crown ratio, the postthinning growth of earlier suppressed trees was retarded clearly.The retarding effect was stronger on spruces than on pines, which supports the results of Vuokila (1970Vuokila ( , 1977)), Hynynen and Kukkola (1989), Eriksson (1990) and Mielikäinen and Valkonen (1991), in which the thinning from above was found to be more convenient for pine stands than for spruce stands.
The aim of this study is to find optimal thinning regimes for Scots pine and Norway spruce mixtures using soil expectation value at 3 % interest rate as the management objective.The study describes a simulation model which uses the distance-dependent single-tree models for diameter and height growth, crown ratio and tree mortality developed by Vettenranta (1999).The models, which take into accout the effect of crown ratio on tree growth and mortality, have not been earlier used in optimizing the management of coniferous mixtures.The optimization problems were solved by using the non-linear method of Hooke and Jeeves (1961).The method for generating coniferous mixtures is described.Six initial stands with different species composition were used in the computations to compare the effect of the structure of the initial stand on the thinning regime and the economic return.Finally, the optimization results are compared with the earlier ones to detect the effect of the new stand simulator on stand management, especially on thinning method.

Initial Stands
The data measured in 9 study plots located in Pieksämäki, Finland (62°20' N, 27°10' E, 110 m a.s.l.), were used for estimating the models for generating the initial stands for computations in this study.Study plots consisted of 133 planted pines and 178 planted spruces.The species composition varied from pure pine to pure spruce stand.The study plots were measured at the age of 15 years for diameter at breast height, height and location of trees (Table 1).The number of stems varied between 1500-2200 ha -1 .
In the simulations, generated stands were used to eliminate the effect of a site quality and the undesirable variation in the structure of initial stands on the stand development and the objective function.The stocking of the generated initial stands was set to be 2000 trees ha -1 , and the proportion of pines varied between 25-85 % of the number of the stems (Table 2).The stands were generated as follows: trees were placed according to the regular distribution with slight random variation and the species of each tree was determined randomly according to the proportion of species.The heights and the diameters were generated by the regression models, which were fitted to the data from study plots (Table 3).For generating the heights of pines the only predictor was the constant (mean height of pines).The height for spruces was predicted by the proportion of pines within a 3-metre radius circle.The tree height was used to predict the tree diameter for both pines and spruces.The residual variation of the models was added randomly to the predictions.The stands were expected to be located on the Myrtillus site type and 6850 km from the equator.

Simulation Model
The development of the initial stands was simulated without stochasticity by the models introduced by Vettenranta (1998;eq. 5-12).The simulation model contained separate diameter and height growth models, as well crown ratio and mortality models for each species.The growth models predict the future 5-year growth using tree dimensions, competition between trees, stand characteristics and location of the stand as predictors.The static crown ratio models use the same predictors.The mortality models remove trees under heavy competition using competition indices and tree heights as predictors.As the predicted variables were used also as predictors in other models, the bias caused by the nonlinear modifications of variables was corrected with the Monte Carlo simulation, in which the independent normally distributed random effect were added several times to the predicted variables (Kangas 1996(Kangas , 1997)).The final result was the mean value of these predictions after the inverse transformation.
The 5-year growth period was simulated as follows: first, the crown ratios were generated with the crown ratio models.Second, the trees under heavy competition were removed by the mortality model.Then, the diameter and height growths were calculated, and after that the stand characteristics with the new tree dimensions.The new crown ratios were calculated with the new tree dimensions and the stand characteristics.
The lenght of the crowns was determined with the crown ratio.If the earlier lower edge of the crown was higher than the new one, the crown ratio was replaced with the earlier one.Finally, the volume of trees and the share of timber assortments were calculated using the functions of Laasasenaho (1982).When competition indices were computed, a temporary buffer zone was generated around the stand plot assuming that the plot was surroundedby similar plots on all sides.
Thinnings were simulated by setting the number of thinnings and the certain basal area limit for each thinning.If this limit was exceeded during the present 5-year growing period, the time step was shortened so that the basal area limit was attained.The selection of trees removed in thinnings is described in detail in the study of Pukkala and Miina (1998).In the thinning algorithm, the trees of each species were divided into three diameter classes (small, medium and large trees) by the number of the stems.The different types of thinnings were implemented by the removal percentages for each diameter classes and the parameters TF(1) and TF(2).The removal was made from the class in which the ratio of the current basal area to the residual basal area was the highest.The parameters TF(1) and TF(2) conduct the order of the tree selection within the specified class.TF(1) denotes the factor which divides the total growing area into the computational growing areas of individual trees according to their diameter.If TF(1) is positive, larger trees receive larger area and conversely.TF(2) denotes the factor which divides the overlapping growing area for trees according to their heights.If TF(2) is positive, a higher tree receives a greater part of the overlapping area and vica versa.If the factors are zero, each tree receives the same computational growing area and the overlapping areas are devided equally for the trees in question.The selection of the trees to be removed starts from the tree which has the highest proportion of the computational growing area to be overlapped.Trees are removed until the removal percentages from the basal areas are reached.

Objective Function
The economic criterion was chosen as a decision criterion when optimal treatment programmes were determined.The investment problem in stand regeneration was formulated using Faustmann´s (1849) theory of land interest (Eq.1).There are several regeneration and treatment alternatives and the one which maximizes the net present value calculated to infinity should be found. max where SEV is the soil expectation value (FIM ha -1 ), u is the regeneration age of a present stand (a), T is time or stand age (a), R is the income from cutting (FIM ha -1 ), i is the interest rate (%/ 100) and C is the expediture (FIM ha -1 ).
When SEV is used as the decision criterion in this form, it is assumed that rotations with equal treatments are repeated to infinity.This is appropriate in this case as all the treatments were simulated from the planting to the final cutting, and the optimal initialmixture for planting was searched for.The weakness of this method is the uncertainty concerning the interest rate and timber prices in the long run.
When the soil expectation value of treatment alternatives were calculated, the minimum diameters and the unit prices of timber assortments as well as the logging costs must be specified.Proportions of timber assortments were calculated according to taper curve and volume functions by Laasasenaho (1982) where minimum diameters were as follows (over bark, cm): The roadside unit prices were used to be able to include the logging cost into the calculations.Waste wood and dead trees were assumed to have no economic value.Other unit prices (FIM m -3 ) were mean prices in Finland in 1996 (Metsätilastollinen vuosikirja 1997): Logging cost depended on the logged volume and the mean size of logged trees (Valsta 1992): ln(c t ) = 5.410 -0.0522 ln(v mean ) + 0.0243 ln(v mean ) 2 -0.445 ln(V r ) + 0.0397 ln(V r ) 2 ln(c c ) = 5.230 -0.0598 ln(v mean ) + 0.0208 ln(v mean ) 2 -0.3840 ln(V r ) + 0.0327 ln(V r ) 2 where c t and c c are the logging costs in a thinning and in a clear felling (FIM m -3 ), respectively, v mean is the mean volume of logged stems (m 3 ) and V r is the total logged volume (m 3 ha -1 ).
The site preparation was expected to cost 750 FIM ha -1 , the planting 1300 FIM ha -1 and the price of plants was 1300 FIM ha -1 (2000 plants per ha).The planting of pine and spruce seedlings was assumed to have the same costs.The future tending of the stand and the supervision of the work were assumed to be 1650 FIM ha -1 as a net present cost (Metsätilastollinen vuosikirja 1997).Thus, the final costs were expected to be 5000 FIM ha -1 at the moment of establishing the stands.

Optimization Method
The optimal thinning regimes were solved by the direct search method of Hooke and Jeeves (1961).The objective function used was not continuous, and therefore any partial derivates based optimization algorithms cannot be used.This does not prevent the utilization of the direct search method, which uses a form of coordinate optimization to solve the problem (e.g.Bazaraa and Shetty 1979).The computer program developed by Osyczka (1984) was modified to solve the optimization problems.The same method and spatial single-tree models have been used earlier to find the optimal treatment scedules for conifer mixtures by Pukkala et al. (1998) and for pure stands by Miina (1996), Pukkala andMiina (1998), andRautiainen et al. (1998).
It is typical for the optimization of stand management that near the optima the solution sur-face is quite flat with numerous minor peaks (Miina 1996, Pukkala andMiina 1998).These peaks represent the treatment decisions concerning individual trees.Therefore, with the nonconvex objective function, it is not certain that the global optimum can be found.Thus, 10 direct search runs were executed using different starting points for the decision variables.The direct search method is sensitive to the initial values of the decision variables.Therefore, for each direct search run, the best of 400 random search points was selected to be the initial point.From a set of 10 local optima, the best one was selected as the solution to the given problem.
The number of decision variables depends on the number of thinnings.The decision variables used in this study were the stand basal area before thinnings, the removal percentages for each diameter class in each thinning, the tree selection parameters for each thinning, and the stand basal area before the final cutting.Thus, the number of decision variables was 9NT + 1, in which NT is the number of thinnings.In random search, the feasible minimum and maximum values for the decision variables were used.The stand basal area ranged between 10-45 m 2 ha -1 , the removal percentages between 0-100 %, and the tree selection parameters from -5 to 5. In the beginning of the direct search, a coordinate search was done with the steps of 0.1 times the range of the decision variables.The search step was reduced if a better solution was not found.The direct search algorithm was stopped when the change steps of the decision variables were less than 0.005 times the range of the decision variables.

Results
The greatest soil expectation value at the 3 % interest rate (8900 FIM ha -1 ) was attained with the initial stand, in which the proportion of pines was 65 % of the number of the stems (Table 4).The regime included two thinnings.The first thinning was conducted at the age of 39 years when the stand basal area was 37 m 3 ha -1 and the dominant height was about 15 m.The post-thinning basal area of 27 m 3 ha -1 was quite high.Pines were thinned from above, so that a half of the largest pines were removed as well as one third of the smallest pines (Table 5).Spruces were thinned from below, so that half of the smallest spruces were removed (Fig. 1).The second thinning was done 8 years after the first thinning, and it was markedly heavier.All the large pines were removed and more than half of the smallest ones.About half of the smallest spruces, some medium-sized spruces and about one third of the largest ones were removed, too.
All the optimal two-thinning regimes of the initial stands had the following features in common: the high stand basal areas before the first thinning (36-40 m 2 ha -1 ) and the high removal percentages of large pines (56-100 %), which decreased the proportion of pines during the rotation.However, the removal percentages of large spruces varied between 0-56 % in the first thinning and between 0-75 % in the second one.There was no pattern for the tree selection parameters.The thinning type is determined mainly by the removal percentages for the diameter classes because the tree selection parameters affect only the order in which trees are removed within the diameter classes.In the optimal one-thinning regimes, the prethinning basal areas were slightly greater (38-41 m 2 ha -1 ) and the post-thinning basal areas (15-18 m 2 ha -1 ) were markedly smaller than in the optimal two-thinning regimes (Table 4).Also, the stand basal area in the final cutting was 1-9 m 2 ha -1 greater and the rotation was 4-7 years shorter than in the two-thinning regimes.The optimal SEVs of the two-thinning regimes were about 5 % higher than those of the one-thinning regimes.Both pines and spruces were removed from all size classes so that the proportion of pines decreased.All large pines were removed in every initial stand (Table 5).
The effect of the density of the initial stand on the SEV was studied by changing the density of the most profitable stand (D) 10 % smaller and greater.At the same time, planting costs were adjusted by 300 FIM ha -1 , respectively.The stand with the lower density (1800 ha -1 ) resulted in almost the same SEV (8840 FIM ha -1 ) as the best outcome, whereas the optimal thinning regime for the stand with the higher density (2200 ha -1 ) was about 5 % less profitable (8430 FIM ha -1 ).
The optimal two-thinning regime was compared to the thinning regime of Scots pine growing on Myrtillus site type in southern Finland presented by Forest Centre Tapio (Luonnonläheinen metsänhoito... 1994).According to the management regime of Forest Centre Tapio, three thinnings were carried out in which the pre-thinning basal areas were 27, 28 and 29 m 2 ha -1 and the post-thinning basal areas were 18, 19 and 20 m 2 ha -1 , respectively.The same removal percentages 95, 45 and 0 in thinnings from below and 30, 10 and 50 in thinnings from above were used for both species.The final cutting was conducted when the stand basal area reached 30 m 2 ha -1 i.e., at the age of 70 years when the stand was thinned from below and at the age of 65 years when the stand was thinned above.With thinnings from below, the loss of SEV was about 30 % (6070 FIM ha -1 ), and with thinnings from above, the loss was about 20 % (7250 FIM ha -1 ).

Discussion
The study presented the simulation-optimization system, which was used to find the optimal regimes for six coniferous mixtures having different initial species composition when the soil expectation value at a 3 % interest rate was maximized.The spatial growth, crown ratio, mortality and thinning models were found to be flexible tools in the simulation of the complex growth dynamics of the mixed forest.The crown ratio models were used to predict the thinning reaction more accurately, so that the development of a stand thinned from above would be more reliable.The Monte Carlo simulation was used to decrease the bias caused by the transformations of dependent and independent variables.However, the results obtained by optimizing a nonconvex objective function using nonlinear programming cannot be verified to be the global optima.
The greatest SEV was attained with the initial stand, in which two thirds of the number of the stems were pines.The optimal thinning regime partly followed the results obtained by Pukkala et al. (1998b).In the present study, the stand basal area raised to quite a high level before the first thinning (37 m 2 ha -1 in the present study and 41 m 2 ha -1 in the study by Pukkala et al. (1998b)).Pines were removed from above and below, but spruces were removed slightly from below only.The second thinning followed the first one quite soon (after 8 years).In the present study, pines were heavely thinned from above and spruces from all size classes, whereas in the comparison study pines were removed more from below and spruces clearly from above.The proportion of pines at the end of the rotation were markedly greater (45 %) in the present study than in the study by Pukkala et al. (1998b).The lower stand basal area before the first thinning in the present study is partly due to the different mortality models and the effect of the crown ratio, which decreases the post-thinning growth of trees the more the lower the crown ratio of trees.The higher unit prices, the faster early development and the optimal structure of the planted stand in the present study may be the reason for the great difference between the soil expectation values of different studies (8900 FIMha -1 in the present study and 6190 FIMha -1 in the comparison study).
The variation among the optimal thinning regimes of the different initial stands is so great that detailed instructions for the treatment of coniferous mixtures cannot be given.Although, some general quidelines can be observed.The early development of pines is markedly faster than that of spruces and the growth rate of spruces is more vigorous than that of pines in the mature forest (Koivisto 1959).Therefore, in the initial stand it is profitable to have quite a great number of pines and then decrease their proportion during the rotation.Early and high net incomes increase SEV the more the higher the interest rate.Thus, it is profitable to grow the stand until several saw-timber pines can be removed in the first thinning, even if it would result in slight openings in the stand.Small and medium-sized spruces and small pines are removed to improve the spatial arrangement, and to avoid the mortality of suppressed trees (Fig. 2).
In the second thinning, all the largest pines and some large spruces with the small relative value growth are removed without caring about the small openings.Small and medium-sized trees of each species are also removed to reach the regular spatial arrangement and a small variation in the size distribution (Fig. 2).The removal of the large trees result in quite a large number of stems (600-800 ha -1 ), even if the stand basal area is not particularly high (< 20 m 2 ha -1 ).
By including the crown ratio into the growth models the reliability of growth predictions after thinning from above was expected to increase.However, the data is mainly consisted of stands thinned from below (common practise in Finnish forests).Although, the data includes small single trees which have been released from competition and the (d or h)/Tg ratio and the crown ratio describe the vitality of those trees, the behaviour of the earlier suppressed trees after the thinning from above is still uncertain.Thus, empirical studies about the thinning reaction after thinning from above is needed.
The simulations made with the stand simulator used in this study indicated that pines promote the height growth of spruces (Vettenranta 1999).The height development of spruces corresponded to that on the Oxalis-Myrtillus site type, while the height development of pines corresponded to that on the poor Myrtillus site type (Koivisto 1959).Due to the fact that the growth of spruce on fertile sites is more rapid than that of pine, and conversely on poor sites, the optimal treatment regime and structure of an initial stand change according to the site quality, as well as according to the geographical location of a stand.Therefore, the optimal treatment schedule cannot be determined in general, but several reasons force to solve the most profitable treatment decisions individually in each case.
The treatment schedules with one-thinning were about 5 % less profitable than those with two thinnings.Thinning regimes with only onethinning can be assumed to be more profitable when using an interest rate higher than 3 % and less profitable when using lower interest rates.In some regimes, about 60 % of the stand basal area were removed which increases the risk of windfall and snow damage.Also thinning from above in a very dense stand may cause severe harvesting damages which were not taken into account in the simulations.However, the results with one-thinning indicate that the soil expectation value is not specially sensitive to the postthinning basal area.The decreased post-thinning growth is compensated by the higher and earlier harvesting income.If the forest owner for some reasons wants to use a higher interest rate than 3 %, then one heavy thinning could be a reasonable alternative.
However, the thinning types applied in the optimal regimes are against the common practise (MMMn päätös 1997) in which thinnings from below are favoured.Also the stand basal areas before the first thinnings are markedly higher and the post-thinning basal areas in mature stands are at lower level than those suggested in the instructions of Forest Centre Tapio (Luonnonläheinen metsänhoito ... 1994).It seems as if the economic goals of forest management conflict with the instructions, and the differences are emphasized if the thinnings are made from above and the interest is taken into account.Of course, the structure and the growth dynamics of mixed forests differ from that of pure stands and that may partly cause the conflict with the instructions which are prepared for the maximum yield of pure stands.
The differences between mixed and pure stands are not necessarily profitable for mixed stands.Pines with a faster initial growth may suffer from quality problems when the canopy layer is not as dense as in a pure stand.On the other hand, the higher density on the lower canopy layer may improve the quality of pine.The optimal tree selection in thinnings may also be very complex, specially when species are not spatially evenly distributed.This results in difficulties in forest management when treatment decisions can no longer be simple general decisions.However, the development of computers and models gives new possibilities for more detailed planning that allows to reach the goals of forest management of individual forest owners better.

Fig. 1 .
Fig. 1.Development of the stand basal area in the optimal one-thinning (A) and two-thinning (B) regimes.

Fig. 2 .
Fig. 2. Stem maps related to the optimal two-thinning regime of stand D at the time of the first thinning (A), second thinning (B) and before the final cut (C).P and S denote pine and spruce, respectively; + and -indicate removed pine and spruce, respectively.

Table 1 .
Main characteristics of the trees in the data used in the stand generation.Basal area weighted mean value in the brackets.

Table 3 .
Parameter estimates and standard errors of the models used in the stand generation.SP-% is the proportion of pines of the number of the stems within 3-metre radius circle centered the subject tree; d is tree diameter (mm) and h is tree height (cm).

Table 4 .
Optimal treatment programmes and soil expectation values (SEV, FIM ha -1 ) at 3 % interest rate for the generated stands.Year is the time since regeneration (a); G and G* are the pre-thinning and post-thinning basal area (m 2 ha -1 ), respectively; SP-% is the proportion of pine of the stand basal area before and after the thinning.

Table 5 .
Optimal removal percentages for small, medium-sized and large pines and spruces, and tree selection parameters.SP-% is the proportion of pines of the number of the stems in the initial stand.Th is the first or the second thinning.