Effects of Catastrophic Insect Outbreaks on the Harvesting Solutions of Dahurian Larch Plantations

Optimal harvesting under pest outbreak risk was studied on a set of even-aged Dahurian larch (Larix gmelinii) stands in northeastern InnerMongolia, China.The effects of catastrophic pest outbreaks caused by the Siberianmoth (Dendrolimus sibiricus) on the economic harvesting plan are compared through both deterministic and stochastic cases. Stand simulation is based on an individual-tree growth system. A scenario approach is applied when simulating the effects of catastrophic pest outbreaks. Insect damage is assumed to be a Poisson process with an average rate of 0.1 per year. One hundred scenarios of insect damage are created using the Poisson process to simulate the distribution of bare land value of each of the optimal regimes. Numerical results show that the optimal rotation is shortened with an increasing probability of a catastrophe.The average bare land values in the stochastic case are approximately 14.8% to 22.9% lower. Numbers of thinnings are decreased for most plots when seeking a highest bare land value, compared to the deterministic optima. If given a constant thinning rate, increasing risk-taking shortens the optimum rotation, as the model set used.


Introduction
Dahurian larch (Larix gmelinii) is one of the main timber species in China.It forms large forests in northeast China, including the northeastern part of the Inner Mongolia autonomous region and the Daxing'an Mountain range.However, because of tree species monocultures and lack of biodiversity, this area has become vulnerable to insect invasions.Disasters caused by the Siberian moth (Dendrolimus sibiricus), one of the most destructive defoliators of larch in the area, have destroyed tens of thousands of hectares (ha) of larch forests and caused immensurable loss to the forest resources in northern China during the past few years [1].These catastrophic insect outbreak disturbances can be viewed as stochastic or random events as we cannot precisely determine if, when, or where they will occur [2][3][4].Earlier studies revealed effects of risk (e.g., fire) on stand yields and developed theoretical models of timber harvesting of forest stands subject to risk (see, e.g., [5,6]).Usually, as stochastic and deterministic processes have different effects on forest ecosystem development and natural disturbances and other uncertainty sources can largely influence forests (see, e.g., [6][7][8][9][10][11][12]), questions have been raised concerning the effects that insect outbreaks will have on forest management solutions.
In our study, a stand management analysis system was designed to support the management of Dahurian larchstands and analyze the influence of catastrophic insect outbreak on harvesting solutions in Aershan area of Inner Mongolia, China.The system consists of stand growth and yield simulator based on individual-tree growth and mortality models and an optimization algorithm which is modified from the study by Hyytiäinen et al. [13] to find the optimal management schedule for a given objective function.The effect of catastrophic insect outbreak was analyzed using the scenario approach where stochasticity is represented by a set of scenarios of the uncontrollable random variables.Our hypothesis is that the probability of insect outbreak affects the financially viable management of Dahurian larch stands.

Study Area and Initial Stands.
The study site is located in Aershan area (47 ∘ 07  -47 ∘ 55  N, 119 ∘ 51  -120 ∘ 57  E, mean elevation 1100 m) of northeastern Inner Mongolia, China.The area belongs to the cold and wet temperate zone with an annual average temperature of −3.2 ∘ C and an average precipitation of 452.11 mm.The total area is 4.83 × 10 5 ha, with forest coverage accounting for 3.45 × 10 5 ha.Forest coverage rate is 71.4% [1].The main forest species in the area are pure Dahurian larch (Larix gmelinii) and white birch (Betula platyphylla) plantations, while Dahurian larch is the only tree species largely planted in the region [1,14].
Measurements from nine inventory plots [1] were used as starting points for stand simulations in our current study.The plots were pure even-aged stands of Dahurian larch.All the plots were treated similarly with human activities, and each had approximately 2000 trees ha −1 at the beginning indicating that a precommercial thinning had occurred at each stand (as was assumed in the bare land value calculation).Biological data, including diameter at breast height (DBH), height, crown width, height of the first live branch, species, and degree of damage (healthy, damaged, dead, and fallen tree), were recorded during the field work.

Growth Simulation.
Stand development simulation is based on individual-tree growth model, mortality model, and ingrowth model [15].The simulation model accepts two alternative ways of describing the stand: either by tree list or by diameter distribution of the representative plot.In the latter case, the system requires the DBH, height, age, and number of trees per hectare of every diameter class, as well as the site index of the stand.
As most of the simulators are based on deterministic relationships, stochasticity may have to be introduced to the models afterwards [16].In our present study, the starting point of the simulation is a deterministic, individual-tree growth model developed by Monserud and Sterba [17] and adjusted for Dahurian larch by Jiang [18].The biological component of the simulator consists of models for (1) basal area increment, (2) crown ratio, (3) height growth, and (4) mortality rate.The models, which are described in detail in the Appendix, represent a common modelling method used in many individual-tree simulators.

Insect Outbreak
Frequency.The effects of epidemic pest outbreak that cause considerable mortality are studied and simulated through the scenario approach where catastrophic insect outbreaks are modelled as random events that damage a part of the growing stock.Stochasticity is represented by a set of scenarios.The random process of insect outbreak is modeled using a Poisson process with outbreak events occurring independently of one another and randomly in time at an average rate of 1/ per unit of time.The insect outbreak probability () is then expressed by where  is the time span between two successive insect outbreaks and  is the average interval between outbreaks.
Previous records [1,19] reveal that outbreak at epidemic level occurred in 1988, 1997, and 2008 separately, which means it happens approximately every 10th year.Figure 1 shows the graph of the exponential distribution for an average interval between outbreaks at epidemic level  = 10 years.
The graph shows that the probability of a pest outbreak occurring increases at a decreasing rate as  increases.The probability that an outbreak will occur within an interval of five years is approximately 0.39.Within a 10-year interval it is approximately 0.63.It is nearly certain that an outbreak at epidemic level will occur within 45 years.
The stochasticity effect was modelled as random events that damage a part of the growing stock.Values of stochastic variables (e.g., timing and intensity of catastrophe in this study) were generated randomly through computer program using Poisson process.We assume that the damaged trees must be harvested by thinning, with a 25% reduction in their stumpage volume and a doubling of the logging costs.These economic parameters were set subjectively, as no data were available from applicable silvicultural conditions.

Economic Data.
The economic data consist of regeneration costs, roadside prices of larch sawlogs and pulpwood, and harvesting costs.A 3.5% interest rate is applied to calculate optimal solutions.The regeneration cost is assumed to be $606.4ha −1 (including land preparation costs, planting and seedling costs, silvicultural operations costs, and tending costs).The roadside prices of larch are $150.5 m −3 and $103.0 m −3 for sawlog and pulpwood, respectively.Merchantable volumes are predicted and computed using the segmented stem taper and volume equations modified from the study by Fang et al. [20].These models are described in the Appendix.
Harvesting costs were kept constant in our study.In detail, the trimming costs employed here are $17.4 m −3 (including both direct and indirect costs) for thinning or selective harvest and $23.6 m −3 (including both direct and indirect costs) for the final harvest.The transportation cost is $0.08 m −3 km −1 , while the average transporting distance to a timber yard is 57 km.As the other fixed costs were partly calculated within the felling costs, the remaining part mainly consists of the loading and storage cost, which is $1.58 m −3 [21].

Objective Function for Optimal Solution.
The optimization model for even-aged larch stands applied here is modified from the study by Hyytiäinen et al. [13].The optimization model used in our study is adjusted so that the effect of a catastrophe damaging a part of the growing stock is incorporated into the model: subject to where  is the bare land value,  = 1, . . .,  harvests,  denotes the timber assortments (for sawlog,  = 1; for pulpwood,  = 2),  is the final harvest,  0 denotes a matrix describing the initial state of the stand,   is stand age at the th harvest,    denotes a matrix describing the stand state before the th harvest at stand age   , ℎ  denotes -dimensional vector as the ratio of trees removed in the th thinning,   is the roadside price of timber assortments (dollars per cubic metre),   denotes the harvested volumes of timber assortments (m 3 ha −1 ),   is volumes of forest lost to the catastrophe (m 3 ha −1 ),    is a random number between 0 and 1, which gives the proportion of trees destroyed by the catastrophe,  is the salvage rate,   is the harvest cost of the th harvest,    is the cost of cleaning up and reforesting the insect-attacked stand area,  0 is establishment and silvicultural costs,  is interest rate, and    denotes the set of admissible thinnings at stand age   .
Gross harvest revenues of the th harvest are summed over ntrees,  timber assortments as products of prices   (dollars per cubic metre), harvested volumes   (m 3 ha −1 ), and salvaged volumes  ⋅   (m 3 ha −1 ).The net harvest revenues are attained by subtracting harvest costs,   and    , consisting of felling and on-site transport costs.The net harvest revenues are discounted to the beginning of the rotation period by a factor (1 + ) −  , where ris the rate of interest.
The bare land value in ( 2) is maximized subject to the stand dynamics that define the stand state just prior to the harvest at age   and the catastrophe scenarios between present (  ) and previous ( −1 ) harvest.This is given as a function of stand state prior to the previous harvest (  −1 ), the intensity of the previous harvest (ℎ −1 ), the possibility of a catastrophe happening and its damage intensity (  −1 ), and the time difference between present   and  −1 .In our study, stand structure is described with  trees ( = 3); each tree is characterized by variables reflecting its current dimensions (age, diameter, height, etc.) and an expansion factor representing the number of trees per hectare.All remaining trees are removed in the final (th) harvest.
The thinning rate equation (5) defines thinning rates for each diameter class.The thinning rate is a piecewise linear function of tree diameter in our study, relative to the smallest and largest stand diameters.The thinning parameters define the thinning rates at the corner points of the piecewise linear function.The thinning definition can be adjusted to simulate different thinning types by changing only a few parameters per thinning (for more details, see [16]).

Method for Finding the Optimal Solution.
Optimizing the management schedule means finding the optimal values for a set of decision variables (DVs), that is, the optimal DV values for each thinning and final harvest [22][23][24][25].Because the number of thinnings is not a continuous variable, schedules with different numbers of thinnings are to be treated as separate optimization problems.The management regime was specified by the number of thinnings and the DVs, chosen as follows.
For each thinning: (a) years since previous thinning (which were set as a multiple of five years from the previous thinning), (b) remaining basal area, and (c) harvest percentage.
For final cuttings: (a) years from the last thinning to the regenerative cut.Noticeably, the criteria for grouping trees in our study are based on their DBH.Tree volumes and stumpage values are thereby computed by diameter classes classified into small stems (with diameters ranging from 4 cm to 14 cm), medium stems (with diameters ranging from 14 cm to 21 cm), and large stems (with diameters greater than 21 cm).
To determine the optimal management schedule, an initial solution (a supposition of the optimal DV combination) was first fed into the computation.By doing this, the value of the objective function with these DV values and the given set of parameters (timber prices, unit costs, and discounting rate) was calculated.The vector of DVs consists of times between forest operations (adjusted by 1 year) and information on how the operations are executed, for example, thinning percentages at different tree diameters (adjusted by 0.1%).Based on the feedback from the simulation system, changes in the DV values were made in an effort to improve the schedule, and the objective function was once again calculated using the new DV values.The optimal management schedule for a given number of thinnings was eventually found, after this search process was repeated several times once a defined convergence criterion was met.The optimization was carried out for 0, 1, 2, 3, 4, and 5 thinnings separately, in order to find the number of thinnings that maximized the bare land value.In order to gain a better economic return, thinning intensity limits for small and medium trees were additionally set at 40%, which a previous study found to be effective for larch plantations in northeastern China [26].

Deterministic Optimal Solutions at 3.5% Interest Rate.
Table 1 shows the deterministic optimal solutions at a 3.5% interest rate for different plots.The symbols in the tables are as follows: DBH is diameter at breast height (cm), S-trees (%) denotes the portion of small trees removed, M-trees (%) denotes the portion of medium trees removed, L-trees (%) denotes the portion of large trees removed, Ro.dg denotes mean diameter at the end of rotation (cm), M.A.I. denotes mean annual increment in cubic metres per year, and SI denotes site index, which is the dominant height in meters at the index age.
Computations of sawlog and pulpwood amounts were based on the diameter classification.For computing the amounts of pulpwood in our study, inferior trees such as small stems and dead or damaged trees (natural mortality or being caused by pest outbreak) were used as pulpwood; the computations of sawlog amounts are based on the volumes of medium and large stems.Thinning type is determined with three variables that define the thinning rates for minimum, middle, and maximum diameters.The thinning rates for other diameters are interpolated.
To maximize the objective function according to various sample plots, two to five thinnings were required.Two thinnings were required on Yiershi 48, three thinnings were required on Dulaer 76, Yiershi 53, Yiershi 57 O-2, and Yiershi 57 O-3, four thinnings were required on Dulaer 80-2, Yiershi 58, and Yiershi 57 O-1, and five thinnings were required on Yiershi 57.The sample plots in our study varied in their initial stand age, site fertility, initial stand density, and geographic location.Prior management treatments also possibly varied.However, bare land values in our study were calculated assuming no prior harvest revenues or previous establishment costs.

Sensitivity Analysis of Varying Numbers of Thinnings.
Deciding on the number of thinnings and the rotation length constitutes separate problems with corresponding numbers of variables.The number of thinnings is therefore a parameter given to the simulation-optimization system at the outset.The maximum bare land values for different numbers of thinnings (i.e., 0, 1, 2, 3, 4, and 5 thinnings) on various plots under optimal rotation lengths are shown in Figure 2.Although bare land value increased with the number of thinnings until the highest was reached, most of the economic gain was achieved with just one thinning.Four thinnings were practically as good as the optimal number (i.e., 3 thinnings) for plots such as Dulaer 76 and Yiershi 53.Five thinnings were practically as good as the optimal number (i.e., 4 thinnings) for Yiershi 58 and Yiershi 57 O-1.For other plots, the bare land values change limitedly with the number of thinnings varying by one.

Effects of Catastrophic Insect Outbreak.
To reveal the effects of catastrophic events or pest outbreaks in our study, scenarios with different probabilities and the intensity of occurring catastrophes were applied.We used 100 scenarios to gain a better picture of the result variability, as the application of fewer scenarios may cause instability of the results.
In even-aged management, one key decision variable is rotation age [27].The same method of scenario approach was used to produce good statistical statements regarding the effect of the rotation age.For this purpose, 100 replications of the simulation at different rotation ages were made, holding all parameters constant except for the string of random numbers (e.g., catastrophe timing and intensity).These 100 scenarios of insect damage were then generated randomly by using the Poisson process and were used to simulate the distribution of present value of each of the optimal regimes obtained from the deterministic analysis.The largest and smallest bare land values, as well as the mean bare land value, the standard deviation of the bare land value, and the standard error of the mean, were observed in these 100 replications with different scenarios for each different rotation.These results are gathered and showed in Table 2.
Table 2 shows that, in these 100 replications, the bare land value ranged from a minimum value to a maximum value.
For example, at rotation age 49, the bare land value for plot Dulaer 76 ranged from a minimum of $6749.0ha −1 to a maximum of $12850.8ha −1 .The mean bare land value was $10218.6 ha −1 , with a standard error of $131.2 ha −1 .The 95% confidence interval of the mean bare land value is 10218.6 ± 2 × 131.2 = ($9956.2ha −1 , $10481.0 ha −1 ).This 95% confidence interval does not contain the bare land value obtained through the deterministic simulation, $11989.5 ha −1 .Thus, the mean results from the stochastic simulation are very different from those of the deterministic simulation.Also, the catastrophes have an effect on the numbers of thinnings when seeking the highest bare land value.indicates that the numbers of thinnings for most plots are decreased in the stochastic case.
In Table 3, the number of thinnings was decreased from four to three for Dulaer 80-2, Yiershi 58, and Yiershi 57 O-1; for Yiershi 53 and Yiershi 57 O-3, the number of thinnings was decreased from three to two; for Yiershi 57, the number of thinnings was decreased from five to four.
The optimal solution for each plot in the deterministic analysis was added to perform the comparisons.The optimal bare land values for both deterministic and stochastic cases are in bold characters.
Table 2 indicates that the mean bare land values from the stochastic simulation (with 100 replications of different scenarios) are 14.8% to 22.9% lower than the deterministic simulation, which corresponds to the different plots.
When the stand was subject to the risk of catastrophe, increasing risk reduced the optimum rotation.This is demonstrated in Figure 3, and we took Dulaer 76 as example, where the deterministic case is compared with three levels of risk, presented by annual probabilities of a catastrophic insect outbreak.Damaged trees are logged by thinnings.Given a constant thinning rate, the optimum rotation was shortened from 53 yrs to 47 yrs when the probability increased from 0 to 0.04.Thinning was brought forward in higher risk cases (Table 4).

International Journal of Forestry Research
The other decision criterion examined in this study was conditioned by a specified bare land value.By adjusting rotation length, the forest manager wishes to maximize the probability of attaining a target bare land value.In this formulation, the decision-maker can increase willingness to take risks by raising the bare land value requirement.When thinning rate was constant, increased risk-taking caused the optimum rotation, that is, the one with highest probability, to decrease (Figure 4).For example, the probability that a simulated bare land value is higher than 10000 USD was 0.63 when the rotation was 49 years.Namely, in 63 of the 100 simulations made, the bare land value was higher than 10000 USD.

Discussion
The main purpose of the present study is to reveal the effects of insect outbreaks on the financial performance of Dahurian larch.To do this, a stand level optimization approach is applied.Stochasticity is introduced in the form of scenarios.
Results in the deterministic case indicate that stands with high mean annual increment such as Yiershi 48, Yiershi 57 O-2, and Yiershi 57 O-3 produced the highest bare land values among all the plots.Many factors affect rotation age, for example, planting density and thinning frequency and intensity [13,25].
Clearcutting would take place around the age of 42 to 58 years on various plots.The mean diameters of all the plots at the end of the rotation period in this deterministic case ranged from 15.9 cm to 21.3 cm, which is different to results received in other studies (see, e.g., [28][29][30]).These studies had longer rotation ages and greater mean diameters (usually >20 cm) at the end of rotation period.On the other hand, Sun et al. [26] and Gao et al. [21] reported similar results of mean diameters at the end of rotation in their research to those reported by our study.The heavy cutting intensity of large trees in previous thinnings and the variance of the initial state and site quality might be key reasons for the difference in mean diameters at the end of rotation.The optimal rotation period for plot Yiershi 57 O-3 is much shorter compared to other plots.Usually larch plantation stands must be 45 years or older to meet clearcutting requirements [21].However, some earlier studies reported that the optimal rotation can be shorter than recommended.For example, Sun et al. [26] claimed that the rotation period of larch is 37-39 years with the first thinning timed 25 years.This result was naturally found using a different discount rate, timber price, and site location.
For the stochastic analysis, the risk of catastrophic pest outbreaks shortens the optimum rotation considerably and decreases the numbers of thinnings, respectively.Among all nine plots, rotation period of Yiershi 53 decreased mostly in stochastic case, by eight years.For Yiershi 57 the risk of catastrophic pest outbreak shortened the optimum rotation period by seven years, for Dulaer 80-2 and Yiershi 57 O-1 both by six years.And for most plots, numbers of thinnings were one time less than the deterministic optima.The increasing catastrophic risks affected the bare land values with variability.The bare land values from the stochastic case were approximately 14.8% to 22.9% lower on average than in the deterministic case.Comparing to the deterministic case, bare land value of Yiershi 57 decreased most by 22.9%, while on the other hand, bare land value of Dulaer 76 decreased only by 14.8%.In addition, with simulating a specified probability (%) we could achieve a bare land value nearly as good as that achieved with the deterministic case when, for example, Dulaer 76, 49 yr rotation period is applied in the stochastic case (deterministic case 53 yr rotation).
Our results point out that initial conditions and management plans may strongly influence the results of including stochasticity into stand level optimization.The source of stochasticity was the occurrence of a catastrophic pest outbreak.To determine the solution with confidence and get a better picture of the variability of the results, 100 scenarios were exerted-a smaller number of scenarios application may cause instability and report incorrect solution.
Information about possible outcome variability has additionally been gained from the stochastic simulation.While the law of averages may be relevant for owners with several stands of this type, such as the owner in our study, it may not be relevant for owners with a single woodlot.Any one of the stochastic simulation outcomes is possible for this owner group, and outcome variability will affect property values and influence managerial decisions.
In our study, the scenarios were formed by generating random outcomes of the defined stochastic process (i.e., catastrophic pest outbreak).Another alternative would be to form scenarios subjectively, so that they would be likely realizations of the future.Furthermore, one of the strengths of the present scenario-based stochastic simulation is that the distribution of the objective function values is known and can be analyzed with respect to relevant variables, by scenario approach.The scenario simulations produce huge amounts of information which can be used to evaluate management strategies.
Quantitative probabilistic risk analyses in our study used mathematical definitions that described the relationship between the expected value of the probability of the event occurring at a particular location and intensity and the net value change given that the event has occurred.Calculating the probabilities for stochastic processes and the effects of catastrophic or expected loss in this manner is advantageous for risk analysis related to disturbance processes.The results and approach offer managers several opportunities because they provide information that can be used to assess risks that are difficult to amend through mitigation or in situations where additional mitigation investments have limited effects.
Forest planning and decision-making usually deal with large forest area, instead of a single stand [31,32].The applicability of the presented system depends on the production objectives of the forest.If the objective is to maximize net income or profitability while considering a stochastic situation or risk of catastrophe, the system can be used directly.It is possible to separately optimize the management of each stand or use the system to produce management instructions for different sites and stand densities in similar circumstances.Stand level optimizations also serve comparative analyses on the effects of economic or biological factors on stand management [32].
We believe that the system presented in our study is a useful contribution towards improved decision-making in the management of Dahurian larch stands under risk of catastrophic insect outbreaks.The system enables us to assess the likelihood and severity of economic consequences when exposed to an invasive species.The system is essentially a stand level decision-support tool, which can be used to simulate alternative management regimes and find economic management strategies.The system can also be used in forestlevel optimization in conjunction with linear programming or other optimization procedures.System reliability will obviously depend on the reliability of the simulation system.It is therefore important to refine and enhance it, for example, by more precisely quantifying economic and ecological damage related to the invasive species, incorporating stochastic structure into the simulation system [33] which leads to valid statistical inference, improved estimation efficiency, and more realistic and theoretically sound prediction.It is proposed that individual-tree modelling methodologies need to characterize and include structured stochasticity in future research.

Figure 1 :
Figure 1: Probability that an insect outbreak at epidemic level occurs within  years, given an average time between outbreaks of  = 10 years.

Figure 2 :
Figure 2: Sensitivity of the maximal bare land value to the number of thinnings on various sample plots.

Figure 3 :
Figure 3: The effect of the increasing risk of a catastrophic insect outbreak on the bare land value of plot Dulaer 76 for different rotations.Legend numbers refer to annual probability of disaster.The discount rate is 3.5%.

Table 2 :
Summary statistics of calculating 100 scenarios of the bare land value under pest risk with different rotations (bare land value, unit: $ USD/ha) and the optimal bare land value under deterministic situation.
Max = maximum bare land value, Min = minimum bare land value, Mean = mean bare land value, SD = standard deviation, and SE = standard error of the mean.

Table 3 :
Effect of catastrophic insect outbreak on the optimum thinning numbers.The discount rate is 3.5%.

Table 4 :
Effect of the increasing probability of a catastrophe ( cat ) on the economic thinning schedule of Dulaer 76.The discount rate is 3.5%.
Figure 4: Probability that a computed bare land value (Dulaer 76) is greater than the specified bare land value target, , for different rotations.The three targets are 9000, 10000, and 11000 USD, with 0.02 probability of catastrophic insect outbreak.Interest rate is 3.5%.