Simulating the effects of wind and snow damage on the optimal management of Norwegian spruce forests

School of Forest Sciences, University of Eastern Finland, PO Box 111, 80101 Joensuu, Finland Sustainable Forest Management Unit (UXFS), Department of Agroforestry Engineering, Higher Polytechnic Engineering School, University of Santiago de Compostela, Campus Universitario, 27002 Lugo, Spain Norwegian Institute of Bioeconomy Research, PO Box 115, 1431 Ås, Norway Forest Sciences Centre of Catalonia, Ctra de St. Llorenç de Morunys, km 2, 25280 Solsona, Spain


Introduction
Forests are facing an increasingly uncertain future due to natural disturbances and changes in their regimes (Lindner et al., 2010;Seidl et al., 2017). Increasing uncertainty due to disturbances has been linked to climate change, variation in forests' interactions with disturbances, lack of knowledge on the anthropogenic suppression of disturbances, and the introduction of invasive species (O'Hara and Ramage, 2013). As a result, disturbances, and their ecological impacts on forests, have been put at the forefront in several studies (e.g. González-Olabarria et al., 2008;Pukkala et al., 2016;Seidl et al., 2017).
Snow and wind disturbances are a common cause of damage in forests (Valinger and Fridman, 2011;Mitchell, 2013) that must be addressed in forest planning. Such disturbances not only have significant ecological effects on forest biodiversity, creating new habitats and modifying the forest carbon balance (Bellone et al., 2017;Zubizarreta-Gerendiain et al., 2017), but also represent an economic threat (Gardiner et al., 2013). Snow and wind damage lead to the unplanned harvesting of trees or to a decision to leave damaged trees unharvested, thus changing the original management plans, and reducing economic revenues and the provision of forest products or services (Gadow, 2000). If the aim is to correctly predict the level of forest services or economic revenues, we need tools to integrate risk into the planning of forest management prescriptions, which will become increasingly relevant in the coming decades .
Although the extent to which climate change will affect the occurrence of snow and wind patterns in boreal forests is currently unclear (Mölter et al., 2016), we can nevertheless address the effect of considering these risks on management (Gadow, 2000). Risk can be considered in forest management by either developing prescriptions that enhance the resistance and resilience of forests, or by anticipating the outputs of forestry under alternative risk scenarios, and using those estimations to select management options that could minimize potential losses.
One approach to developing forest management alternatives that include risk is to use a computational representation of the forest system (e.g. Valinger and Fridman, 1999). Forest dynamics are simulated with models created using empirically collected data (Seidl, 2016). Through simulations, we are able to produce detailed representations of forest dynamics, and obtain results that were never measured when the forest dynamics were modelled (Winsberg, 2010). For example, we can investigate the effect of completely new management options on the desired forest functions, such as the effect that changes in thinning intensity or rotation length will have on forest production. When forest dynamics are simulated with one of the processes being stochastic, the simulation provides a range of possible outcomes. As a result, the manager can analyse the risks and uncertainty associated with alternative management options (Pukkala and Kangas, 1996;González-Olabarria et al., 2017;Zubizarreta-Gerendiain et al., 2017).
Finding the optimal management strategy under uncertain circumstances is a complex task, as it requires the consideration of multiple sources of uncertainty. These uncertainties might result from future changes in harvesting costs and prices (Leskinen and Kangas, 1998;Pukkala and Kellomäki, 2012), management objectives or various biotic and abiotic risks (Heinonen et al., 2009). Under such circumstances, the management aim might be to minimize risk (Heinonen et al., 2009), to see the effect of risk considerations on optimal management (e.g. with timber production as the main objective; Zeng et al., 2007), to maximize the expected timber revenues under risk, to increase forest resilience or to create less vulnerable forest structures (Zubizarreta-Gerendiain et al., 2017).
The selection of an optimal management that considers risk is preceded by estimation of the loss caused by potential damage versus the cost of action to prevent the damage (Hanewinkel et al., 2011). These estimations are done using models that predict the probability and degree of damage. In Norway, there have been studies on the optimal management strategies for steady timber production at the regional level (Hofstad, 1991), timber production potential at the country level (Hoen et al., 2001;Eid et al., 2002) and optimal management considering carbon flows (Raymer et al., 2009;Buongiorno et al., 2012). However, none of the previous studies have included the risk of damage due to abiotic agents such as snow and wind in the optimization of stand management, even though wind and snow disturbances are a common cause of damage in Norway (Díaz-Yáñez et al., 2016).
This study developed and demonstrated an approach to integrating snow and wind damage into forestry decision-making. It also analysed the implications of the risk of damage and the effect of uncertain silvicultural costs and timber prices on optimal management.

Stand dynamics simulation
Different models were used to simulate forest dynamics (Table 1), consisting of tree growth, recruitment and mortality. Additional models were employed to calculate tree height, biomass and volume. The risk of wind and snow damage was calculated using predictive damage models ( Figure 1).
The simulation of growth and forest regeneration was based on the models of Bollandsås et al. (2008). These models predict the 5-year diameter growth of a tree and the expected number of recruits. All the recruited new trees were assigned a diameter of 5 cm, as all the ingrowth models were developed for a 5-cm-diameter limit. The simulation was done at tree level, using representative trees; that is, each tree represented a certain number of trees per ha. A representative tree was removed from the simulation when it represented fewer than 0.5 trees per ha.
A limitation of the mortality model of Bollandsås et al. (2008), from the point of view of our study, was the lack of information concerning the cause of mortality. Therefore, it was not possible to tell how much of the predicted mortality was due to competition and how much was due to damage. Because we used separate models for damage, we fitted a new model for competition-induced mortality during five years, using the same predictors as in the original model of Bollandsås et al. (2008) (see supplementary material). The dataset consisted of trees present in plots with no disturbance records, according to the Norwegian National Forest Inventory (NFI) protocols (Landsskogtakseringens, 2008). The data were obtained from the 8th, 9th and 10th NFIs. All the predictors had to be significant at the 0.05 level, and the residuals had to indicate a non-biased model. To assess the performance of the fitted models, the area under the receiver operating characteristic curve was calculated (often referred to as AUC).
The tree heights were estimated for each 5-year period during the simulation by using the ordinary least square models of Sharma and Breidenbach (2014). These models use the tree diameter at breast height (DBH) as the predictor. The individual tree heights were then used to calculate the mean height predictor of the damage models (Díaz-Yáñez et al. 2017). However, the damage models were constructed using the mean height estimates available in the NFI datasets (Landsskogtakseringens, 2008), and these mean heights are different from those based on the models of Sharma and Breidenbach (2014). Therefore, in the simulation of damage, we used a correction factor to calibrate the mean height estimates to the same level as used in the damage models.
The tree volumes were estimated using existing models: Vestjordet (1967) for spruce, Braastad (1967) for pine and Braastad (1966) for birch and other broadleaves. These models predict the stem volume under bark, based on DBH and tree height.
The initial stand conditions for the simulation (Table 2) were based on a randomly selected plot from the NFI data (8th, 9th and 10th NFIs) representing spruce-dominated forest. The pool from which this initial plot was selected was restricted to: 1) well-stocked young stands with more than 1500 trees ha −1 , corresponding to the youngest development classes, according to the NFI protocols (Landsskogtakseringens, 2008); 2) productive forest land; 3) plots with no silvicultural treatments during the past 5 years; and 4) location within the 25 per cent and 75 per cent thresholds (expressed as quartiles from the distribution of young spruce stands) of slope, altitude and latitude (Table 2).

Risk simulation
The models for predicting snow and wind damage occurrence at the tree (Díaz-Yáñez et al., 2017) and stand (Díaz-Yáñez et al., 2018) levels included both site variables (latitude or altitude) and stand variables (stand density, mean diameter, mean height). The rate of snow and wind damage was predicted from stand basal area, mean tree height, mean diameter and average slenderness of the trees. The models were different for uprooted and broken trees. The damage probabilities obtained from the models of Díaz-Yáñez et al. (2018) were converted to 5-year probabilities. The predicted damage was assigned to different diameters following the observed diameter distributions of damaged trees by dominant stand species (see supplementary material).
Damage was simulated using either deterministic or stochastic approach ( Figure 2). In the deterministic approach, the predicted probability of the stand to be damaged (damage occurrence) was multiplied by the predicted probability of any tree in the stand to be uprooted and broken (damage rate). The resulting value was multiplied by the number of remaining trees in the stand to obtain the number of damaged trees. In the stochastic approach, randomness was simulated (1) at plot level or (2) at both plot and tree level. Stochastic damage was simulated by comparing the predicted probability of damage to a random number drawn from a binomial distribution. In addition, we analysed the impact of applying fixed probabilities for damage occurrence at the stand-level Simulating the effects of wind and snow damage (deterministic damage). The 5-year probabilities were 0.2, 0.5 and 0.8, referred to as low, medium and high probability of damage, respectively.
We also calculated the uncertainty associated with each of four optimal management schedules-(1) no damage, (2) deterministic damage, (3) stochastic damage at plot level, and (4) stochastic damage at plot and tree levelby simulating them 6000 times with stochastic damage at the plot and tree levels. The four optimal management schedules were developed with average stand establishment cost (1500 EUR ha −1 ) and a 2 per cent discount rate.

Economic parameters
The timber values were calculated from the estimated pulpwood and sawn wood proportions of the harvested trees. The pulpwood proportion was estimated using the equations developed for spruce and pine (Blingsmo and Veidahl, 1992) and assuming that 80 per cent of the birch and other broadleaf volume went to pulpwood production (and the remaining 20 per cent was left in the forest). The equations from Blingsmo and Veidahl (1992) assumed perfect trees, and therefore wood waste was not considered. We imposed a diameter limit of 13 cm, assuming that smaller trees did not provide any revenue. We also set a maximum limit of 60 cm (capacity limit of sawmills), beyond which all the volume was pulpwood. The average prices per m 3 were obtained from Statistics Norway (2017), collected in 2016 for spruce, pine and broadleaves/birch. The average sawn wood prices were 50.56, 48.77 and, 32.47 EUR m −3 , and the average pulpwood prices were 23.4, 21.4 and, 24.4 EUR m −3 for spruce, pine and birch, respectively.
Harvesting and forwarding costs were calculated using the productivity functions from Dale et al. (1993) and Dale and Stamm (1994) (Table 3), assuming 100 EUR h −1 in forwarding and 156 EUR h −1 in harvesting, following Bergseng et al. (2013). We applied an adjusted time by slope, following Bergseng et al. (2013). The equation for transport productivity (h m −3 ) presented in Dale and Stamm (1994) was based on the forwarding distance. We assumed a distance to road of 692 m (mean forwarding distance; Antón-Fernández and Astrup, 2012) and a volume per load of 15 m 3 . The total costs of harvesting and forwarding were constrained within a reasonable range: in thinning, between 8 and 20 EUR m −3 , and for clear felling, between 5 and 15 EUR m −3 . A fixed cost of 100 EUR per cutting was also included.
The economic outcomes of the management schedules depend on discount rate (opportunity cost), stand establishment costs and timber prices, all of which are uncertain variables. To assess their impact on optimal management, we performed a sensitivity analysis for each of them, where the stand establishment costs were 1000, 1500 and 2000 EUR ha −1 , the discount rates were 1, 2 and 3 per cent and the timber prices were increased or decreased by 20 per cent compared to the baseline prices given above.

Optimization
The objective of the optimization was to maximize the economic revenue using the land expectation value (LEV) depending on decision variables and variables subject to uncertainty. The decision variables, or controlled factors, defined the management schedule. The purpose of optimization was to identify the optimal even-aged management regime, including decisions on the length of the rotation and the thinning regime. The LEV is the present value per unit area (EUR ha −1 ), assuming a sequence of even-aged forest rotations identical to the first rotation.
In order to calculate the LEV, we first estimated the net present value for the first rotation (R 1 ) (Equation 1).  Simulating the effects of wind and snow damage The first term of Equation (1), -C 0 , corresponds to the stand establishment costs. The second term is the sum of the present values of the intermediate costs, C t , and revenues, P it Y it , in year t, where Y it is the yield (the timber volume cut in the thinning) by species and assortment, i, and P it is the timber price of species and assortment, i. The term r is the discount rate. The last term corresponds to the present value of the final clear fell (at rotation time R).
After calculating the net present value (NPV) of the first rotation, the NPV R1 was converted into land expectation value by applying the following formula for infinite periodic payment: Management was optimized by determining the best combination of decision variables: thinning years, rotation length, thinning intensities and thinning type. Thinning intensities and type were calculated using a thinning intensity curve (Pukkala et al., 2014;Pukkala, 2015;Jin et al., 2017). This curve defines the proportion of harvested trees of a certain diameter. The curve has two parameters, which define the type and intensity of thinning. In summary, the optimized variables were: number of five-year intervals between initial age and the first thinning, between subsequent thinnings (maximum of five thinnings) and between the last thinning and final harvest; and the parameters of the thinning intensity curve for each thinning.
For the optimization, we used differential evolution, a populationbased method which has previously shown good results (Arias-Rodil et al., 2015). This technique uses a population of solution vectors that include the values of decision variables. Different solution vectors are combined to create new solution vectors. If a new solution vector is good enough in terms of the LEV, it replaces one of the previous solution vectors. This process is continued for several iterations, and the best solution vector after the final iteration is taken as the optimal solution. A more detailed description of the algorithm is presented in Pukkala (2009). Based on preliminary analyses, we used 50 solution vectors as the population size and 50 iterations. The other parameter values of differential evolution were the ones recommended in Jin et al. (2018).
Several penalties and restrictions were used in the optimization process (Table 3). The penalty functions were strong enough to prevent illogical solutions and extrapolation beyond the ranges of the models, with: • The simulation being run for a maximum of 205 years, but with a penalty function being applied when the simulation period was longer than 100 years. • Thinning intensity being penalized when the proportion of removed basal area was above 0.7. • Thinnings being restricted to stands with mean diameters above 13 cm. • A penalty being applied for rotations shorter than 50 years. • A minimum interval of 5 years between cuttings, being consistent with the projection period of the growth simulator.

Results
The inclusion of snow and wind risk had an effect on the optimal management strategy (Figure 3). Considering risk either shortened or maintained the rotation length, keeping it equal for most of the optimization cases. Increasing discount rates meant a decrease in the stand basal area and a shorter rotation length. The growing stock volume was lower with risk, especially towards the end of the rotation.

Forestry
Among the management schedules that were optimized under risk, the mean volume losses were 0.25 m 3 ha −1 yr −1 , with a standard deviation of 0.29 (mean of the 27 optimal schedules, excluding the options that modified the probability of occurrence and prices). Simulating the optimal management obtained without considering the risk of damage overestimated the LEVs. The maximum overestimation occurred when the stand establishment cost was 2000 EUR ha −1 and the discount rate was 2 per cent; in this case, the differences in the LEV, with and without risk, were 25 per cent.
Increasing damage occurrence meant a reduction in the rotation length when using the baseline parameters (stand establishment costs of 1500 EUR ha −1 and a 2 per cent discount rate) and a deterministic approach to consider risk (Figure 4). The results showed that setting the 5-year probability of snow and wind risk occurrence at 0.2 gave a similar optimal management as obtained when the risk was predicted using the models. On the other hand, when the probability of damage occurrence was set at a higher level, the optimal management changed: at 0.5, one thinning was applied and the rotation length was 75 years, and at 0.8, the damage was equal to the basal area growth in the stand.
Changing timber prices ( Figure 5) and stand establishment costs ( Figure 6) affected the optimal management. Considering risk using the deterministic approach and by reducing prices by 20 per cent resulted in a longer rotation length, and increasing the prices by 20 per cent reduced the rotation length. Higher stand establishment costs delayed investment in the next tree generation, making the rotation longer. Management considering risk tended to have lower stand volumes. This was most visible when the stand establishment costs were 1500 EUR ha −1 .
Each management schedule had a certain degree of uncertainty as a result of stochastic damage. Figure 7 shows four different optimal management prescriptions with their associated uncertainty. For management considering risk with a deterministic approach, the 95 per cent confidence limits of the LEV were 890-1110 EUR ha −1 . The difference between the mean value obtained from considering risk with a deterministic approach and the optimal management without considering damage was 3.2 per cent.

Discussion
Addressing the effects of wind and snow damage in forest management is a complex task that requires appropriate methodological approaches. However, ignoring the effect of damage in management results in overestimated revenues and suboptimal management, as shown by our results and in other studies (e.g. Zubizarreta-Gerendiain et al., 2017).
Our simulations showed that considering the risk of snow and wind damage tends to shorten the optimal rotation period, as suggested by earlier studies, where longer rotations tend to increase damage (e.g. Wallentin and Nilsson, 2014). In those cases where the rotation length was longer when considering damage, the optimal management approach included a greater number of thinnings or greater thinning intensity, which would explain why the rotation length was longer. Our study only considered timber production. If we included other ecosystem services, such as biodiversity and carbon, it is likely that the effect Table 3 Models used for the economic evaluation, with the restrictions or penalties applied to them. of the thinning intensity curve Thinning intensity is not applied to those stands with mean dbh < 13cm.
There is a penalty for thinnings above 70% of the basal area.

None
The valid area of the function is D=130-450 mm and H=140-280 dm.
Trees with diameters above 60 cm do not produce sawn wood. Trees bellow 13 cm diameter do not provide any revenue.
Harvesting and forwarding costs Dale et al. (1993), Dale and Stamm (1994) HarvProd ( Simulating the effects of wind and snow damage of risk on the rotation length would be different (Liski et al., 2001). Management considering the risk of snow and wind damage tends to have lower growing stock volumes towards the end of the rotation; increasing discount rate (higher opportunity costs for standing trees) had a similar effect. When the stand was managed under risk, it was better to remove trees earlier and decrease stand density when the trees were economically valuable. The same effect has previously been observed when optimizing management under the risk of fire (González-Olabarria et al., 2008).
Considering the effect of disturbances with a deterministic approach does not fully describe the uncertainty associated with the management schedule. We need the mean values of the disturbance effects, but also their variation in order to better evaluate the impact of the disturbance (McCarthy and Burgman, 1995). The use of the stochastic approach reflected this variation, which is more informative for adaptive management, and makes it possible to consider the risk attitudes of the forest landowner. Whereas a risk-averse landowner might prefer less variation, the risk-taker would select management options with greater variation that may result in high revenues with low probability, under exceptionally favourable states of nature.
The results showed that the estimated LEV variation was similar for the four optimal schedules, corresponding to the different approaches used to consider risk, when all schedules were simulated many times so that the damage was stochastic. Since the distributions of the LEV were fairly similar for all schedules, it can be concluded that the risk attitude of the decisionmaker will not affect the ranking of these four schedules. But even in those cases where the ranges of variation are similar, it is still important to know the associated uncertainties for a more effective risk management (Blennow et al., 2014). In our results, there was an improvement in the mean LEV when using the optimal management that considered risk in optimization (with a deterministic approach) compared to the management that was optimized without the risk of damage. On the other hand, management based on stochastic optimization provided lower, or very similar, mean LEV values, probably reflecting that Figure 3 Optimal rotation length and basal area dynamics across time for 1 per cent, 2 per cent and 3 per cent discount rate, stand establishment costs of 1500 EUR ha -1 , baseline timber prices and considering the risk of damage using four different approaches: 1) no risk of wind and snow damage (No damage); 2) deterministic simulation of damage (Damage); 3) stochastic simulation of damage at plot and tree level (Stochastic (pt)); and 4) only at plot level (Stochastic (p)).

Figure 4
Effect of the probability of the occurrence of risk damage on optimal management when the stand establishment costs are 1500 EUR ha −1 and the discount rate is 2 per cent. The lines represent the optimal management with a deterministic simulation of risk with three occurrence probabilities: low (0.2), medium (0.5) and high (0.8).

Figure 5
Effect of changes in the timber prices on management when the stand establishment costs are 1500 EUR ha -1 and the discount rate is 2 per cent. The lines represent the optimal management with a deterministic risk simulation.
Forestry the stochastic simulation of risk resulted in a more complicated optimization problem, making it more difficult to determine the true optimal solution.
In this study, we used a simulator specifically created for Norwegian conditions and an optimization approach to analyse alternative management schemes. The developed simulator presented some advantages over earlier simulators (Eid and Hobbelstad, 2000;Gobakken et al., 2008) in terms of methodological approach, as well as the integration of the most recent models for forest growth (Bollandsås et al., 2008;Sharma and Breidenbach, 2014) and forest damage (Díaz-Yáñez et al., 2017. However, the models used could also be considered to be a limitation, since part of the models were taken from rather old studies. For instance, one drawback concerning the estimation of revenues (Blingsmo and Veidahl, 1992) was the use of models that were developed for a timber market situation that may no longer correspond to the current markets. A better alternative would be using current market prices and a taper function for calculating assortment volumes. However, there are no taper models in Norway, and using the functions from other countries may also lead to errors because these models do not necessarily represent Norwegian conditions. Furthermore, timber prices change through time, increasing the uncertainty in management decisions (Leskinen and Kangas, 1998). We took this aspect into consideration by evaluating the effect of changing timber prices on the optimal management. The results showed that a change in timber price had an effect on optimal rotation length.
In addition to the lack of up-to-date models, the combined use of many models made it difficult to evaluate the reliability of the simulation predictions, due to the propagation of errors (Winsberg, 2010). However, correlations of the errors in different models should be known in order to properly simulate variation not explained by the models. These correlations are largely unknown because the models of the simulator were obtained from different studies, and each study reports the residual variation of only one model (i.e. Bollandsås et al., 2008;Díaz-Yáñez et al., 2018). Figure 6 Effect of changes in stand establishment costs on optimal management for a 2 per cent discount rate. The lines represent the optimal management considering the risk of damage with four different approaches (see Figure 3).

Figure 7
Distribution of LEV (land expectation value) when the risk is considered stochastic at plot and tree level for the four optimal management schedules obtained using four different approaches to consider risk (calculated with a 2 per cent discount rate and 1500 EUR ha -1 of stand establishment costs).
Simulating the effects of wind and snow damage Concerning the representativeness of the results, the analysis was based on a single stand. The main aim of the selected initial stand was to be representative of spruce-dominated young stands. It is obvious that the structure, location and species composition of the stand will affect the results of the simulations. For instance, latitude, altitude and slope greatly influence both the growth and damage predictions (Bollandsås et al., 2008;Díaz-Yáñez et al., 2018). A way to deal with this variation, although beyond the scope of this study, would be to produce a landscape-level analysis using a variety of conditions (Heinonen et al., 2011;Senf and Seidl, 2017).
Regarding the approach for finding optimal management, previous studies considering risk have used direct search optimization methods, such as Hooke and Jeeves (González et al., 2005;González-Olabarria et al., 2008;Möykkynen and Pukkala, 2010). Although all the optimization approaches work adequately, it has been proven that, for more complicated cases where many decision variables are involved, it may be worthwhile using population-based methods, especially differential evolution and particle swarm optimization (Arias-Rodil et al., 2015;Jin et al., 2018).
Several clear conclusions can be drawn from this study. Considering risk of damage from snow and wind has management implications, such as shortening rotation lengths, reducing stand density and the value of the growing stock towards the end of the rotation. The current results improve our understanding of stand management under risk, and highlight the importance of considering natural disturbances when making management decisions.

Supplementary data
Supplementary data are available at Forestry online.