1
Effects of Catastrophic Insect Outbreaks on the Harvesting Solutions of Dahurian
2
Larch Plantations
3
4
Qi Jin1,*, Lauri Valsta1, Kari Heliövaara1, Jing Li2,3, Youqing Luo2 and Juan Shi2
5
6
1
Department of Forest Sciences, P.O. Box 27, FI-00014, University of Helsinki, Finland
7
8
2
9
Forestry University, Beijing 100083, P.R.China
Key Laboratory for Silviculture and Conservation of Ministry of Education, Beijing
10
11
3
China Forestry Publishing House, Beijing 100009, P.R. China
12
13
14
*Corresponding author.
E-mail address: [email protected]
15
16
17
18
19
20
21
22
0
23
Abstract:
24
25
Optimal harvesting under pest outbreak risk was studied on a set of even-aged Dahurian larch
26
(Larix gmelinii) stands in northeastern Inner Mongolia, China. The effects of catastrophic
27
pest outbreaks caused by the Siberian moth (Dendrolimus sibiricus) on the economic
28
harvesting plan are compared through both deterministic and stochastic cases. Stand
29
simulation is based on an individual-tree growth system. Combined with the simulation
30
system, the modified optimisation function is able to find the optimal number, intensity and
31
type of thinnings, subject to given rotation lengths. The objective is the bare land value with a
32
3.5% discount rate. Optimal rotation age may vary considerably (42–58 yrs) for different
33
larch stands depending on the initial stand state in the deterministic analysis. A scenario
34
approach is applied when simulating the effects of catastrophic pest outbreaks. Insect damage
35
is assumed to be a Poisson process with an average rate 0.1 per year. One hundred scenarios
36
of insect damage are created using the Poisson process and used to simulate the distribution
37
of bare land value of each of the optimal regimes obtained from the deterministic analysis.
38
Numerical results show that the optimal rotation is shortened with the probability of a
39
catastrophe. The average bare land values in the stochastic case are approximately 14.8% to
40
22.9% lower. Numbers of thinnings are decreased for most plots when seeking a highest bare
41
land value, compared to the deterministic optima.
42
43
Key words: Dahurian larch (Larix gmelinii), Individual-tree growth model, Scenario approach,
44
Risk, Stand-level Optimization, Siberian moth (Dendrolimus sibiricus)
1
45
1. Introduction
46
47
Dahurian larch (Larix gmelinii) is one of the main timber species in China. It forms large
48
forests in northeast China, including the northeastern part of the Inner Mongolia autonomous
49
region and the Daxing'an Mountain range. However, because of tree species monocultures
50
and lack of biodiversity, this area has become vulnerable to insect invasions. Disasters caused
51
by the Siberian moth (Dendrolimus sibiricus), one of the most destructive defoliators of larch
52
in the area, have destroyed tens of thousands of hectares (ha) of larch forests and caused
53
immensurable loss to the forest resources in northern China during the past few years (Li et
54
al., 2009). These catastrophic insect outbreak disturbances can be viewed as stochastic, or
55
random events as we cannot precisely determine if, when or where they will occur (Peltola et
56
al., 1999; Forsell, 2003; Hanewinkel et al., 2010). Earlier studies revealed that effects of risk
57
(e.g. fire) on stand yields and developed theoretical models of timber harvesting of forest
58
stands subject to risk (e.g. Reed, 1984; Reed and Errico, 1985). Usually, as stochastic and
59
deterministic processes have different effects on forest ecosystem development, and natural
60
disturbances and other uncertainty sources can largely influence forests (e.g. Reed and Errico,
61
1985; Wang et al., 2006; Kerns and Ager, 2007; Jenkins et al., 2008; Moore et al., 2008;
62
González-Olabarria and Pukkala, 2011; Mitchell, 2013), questions have been raised
63
concerning the effects that insect outbreaks will have on forest management solutions.
64
65
In our study, a stand management analysis system was designed to support the management
66
of Dahurian larch stands and analyse the influence of catastrophic insect outbreak on
2
67
harvesting solutions in the Aershan area of Inner Mongolia, China. The system consists of a
68
stand growth and yield simulator based on individual-tree growth and mortality models, and
69
an optimisation algorithm which is modified from Hyytiäinen et al. (2005), to find the
70
optimal management schedule for a given objective function. The effect of catastrophic insect
71
outbreak was analysed using the scenario approach where stochasticity is represented by a set
72
of scenarios of the uncontrollable random variables. Our hypothesis is that the probability of
73
insect outbreak affects the financially viable management of Dahurian larch stands.
74
75
2. Materials and methods
76
77
2.1 Study area and initial stands
78
79
The study site is located in the Aershan area (47°07′ –47°55′N, 119°51′–120°57′E, mean
80
elevation 1100 m) of northeastern Inner Mongolia, China. The area belongs to the cold and
81
wet temperate zone with an annual average temperature of −3.2 °C and an average
82
precipitation of 452.11 mm. The total area is 4.83×105 ha, with forest coverage accounting for
83
3.45×105 ha. Forest coverage rate is 71.4% (Li et al., 2009). The main forest species in the
84
area are pure Dahurian larch (Larix gmelinii) and white birch (Betula platyphylla) plantations,
85
while Dahurian larch is the only tree species largely planted in the region (Li et al., 2009;
86
Wang et al., 2009).
87
88
Measurements from nine inventory plots (Li et al., 2009) were used as starting points for
3
89
stand simulations in our current study. The plots were pure even-aged stands of Dahurian
90
larch. All the plots were treated similarly of human activities, and each had approximately
91
2000 trees ha-1 at the beginning indicating that a pre-commercial thinning had occurred at
92
each stand (as was assumed in the bare land value calculation). Biological data, including
93
diameter at breast height (DBH), height, crown width, height of the first live branch, species,
94
and degree of damage (healthy, damaged, dead, fallen tree) were recorded during the field
95
work.
96
97
2.2 Simulation
98
99
2.2.1
Growth simulation
100
101
Stand development simulation is based on individual tree growth model, mortality model and
102
ingrowth model (Huber, 2010). The simulation model accepts two alternative ways of
103
describing the stand: either by tree list or by diameter distribution of the representative plot.
104
In the latter case, the system requires the DBH, height, age and number of trees per hectare of
105
every diameter class, as well as the site index of the stand.
106
107
As most of the simulators are based on deterministic relationships, stochasticity may have to
108
be introduced to the models afterwards (Valsta, 1992a). In our present study, the starting
109
point of the simulation is a deterministic, individual-tree growth model developed by
110
Monserud and Sterba (1996) and adjusted for Dahurian larch by Jiang (2009). The biological
4
111
component of the simulator consists of models for (1) basal area increment, (2) crown ratio,
112
(3) height growth, and (4) mortality rate. The models, which are described in detail in the
113
Appendix, represent a common modelling method used in many individual-tree simulators.
114
115
2.2.2
Insect outbreak frequency
116
117
The effects of epidemic pest outbreak that cause considerable mortality are studied and
118
simulated through the scenario approach where catastrophic insect outbreaks are modelled as
119
random events that damage a part of the growing stock. Stochasticity is represented by a set
120
of scenarios. The random process of insect outbreak is modeled using a Poisson process with
121
outbreak events occurring independently of one another and randomly in time at an average
122
rate of 1/a per unit of time. The insect outbreak probability P(T) is then expressed by:
123
124
P(T) = 1 – exp(– T/a)
(1)
125
126
where T is the time span between two successive insect outbreaks, and a is the average
127
interval between outbreaks. Previous records (Forest Diseases and Insect Pests Control and
128
Quarantine Station, 2006; Li et al., 2009) reveal that outbreak at epidemic level occurred in
129
1988, 1997 and 2008 separately, which means it happens approximately every 10th year.
130
Figure 1 shows the graph of the exponential distribution for an average interval between
131
outbreaks at epidemic level a = 10 years. The graph shows that the probability of a pest
132
outbreak occurring increases at a decreasing rate as T increases. The probability that an
5
133
outbreak will occur within an interval of five years is approximately 0.39. Within a 10-year
134
interval it is approximately 0.63. It is nearly certain that an outbreak at epidemic level will
135
occur within 45 years.
136
137
Figure 1.
138
139
The stochasticity effect was modelled as random events that damage a part of the growing
140
stock. Values of stochastic variables (e.g. timing and intensity of catastrophe in this study)
141
were generated randomly through computer program using Poisson process. We assume that
142
the damaged trees must be harvested by thinning, with a 25% reduction in their stumpage
143
volume and a doubling of the logging costs. These economic parameters were set subjectively,
144
as no data were available from applicable silvicultural conditions.
145
146
2.2.3
Economic data
147
148
The economic data consist of regeneration costs, roadside prices of larch sawlogs and
149
pulpwood and harvesting costs. A 3.5% interest rate is applied to calculate optimal solutions.
150
The regeneration cost is assumed to be $606.4 ha-1 (including land preparation costs, planting
151
and seedling costs, silvicultural operations costs, and tending costs). The roadside prices of
152
larch are $150.5 m-3 and $103.0 m-3 for sawlog and pulpwood, respectively. Merchantable
153
volumes are predicted and computed using the segmented stem taper and volume equations
154
modified from Fang et al. (2000). These models are described in the Appendix.
6
155
156
Harvesting costs were kept constant in our study. In detail, the felling costs employed here are
157
$17.4 m-3 (including both direct and indirect costs) for thinning or selective harvest, and
158
$23.6 m-3 (including both direct and indirect costs) for the final harvest. The transportation
159
cost is $0.08 m-3 km-1, while the average transporting distance to a timber yard is 57 km. As
160
the other fixed costs were partly calculated within the felling costs, the remaining part mainly
161
consists of the loading and storage cost, which is $1.58 m-3 (Gao et al., 2009).
162
163
2.3 Objective function for optimal solution
164
165
The optimisation model for even-aged larch stands applied here is modified from Hyytiäinen
166
et al. (2005). The optimisation model used in our study is adjusted so that the effect of a
167
catastrophe damaging a part of the growing stock is incorporated into the model:
168
max
{hu ,tu ,u=1,…,k; Z0 }
V=
∑ku=1 ��∑ni=1 ∑2j=1 pj ∙ �gj �Ztu , hu � + s ∙ Lj �Ztu , btu �� − Cu �Ztu , hu � − C′u �Ztu , btu �� (1 + r)−tu − C0 (Z0 )
1 − (1 + r)−tk
169
(2)
170
171
172
173
174
subject to
Ztu = f�Ztu−1 , hu−1 , t u − t u−1 , btu � , ∀ u = 1,..., k,
tu ≤ tu+1 ,
hu ∈ σtu ,
∀ u = 1,..., k−1,
∀ u = 1,..., k,
(2.1)
(2.2)
(2.3)
7
175
176
where V is the bare land value, u = 1,…, k harvests, j denotes the timber assortments (for
177
sawlog, j = 1; for pulpwood, j= 2), k is the final harvest, Z0 denotes a matrix describing the
178
initial state of the stand, tu is stand age at the uth harvest, Ztu denotes a matrix describing
179
180
the stand state before the uth harvest at stand age tu, hu denotes a n-dimensional vector as the
ratio of trees removed in the uth thinning, pj is the roadside price of timber assortments
182
(dollars per cubic metre), g j denotes the harvested volumes of timber assortments (m3 ha−1 ),
183
and 1, which gives the proportion of trees destroyed by the catastrophe, s is the salvage rate,
181
Lj is volumes of forest lost to the catastrophe (m3 ha−1 ), btu is a random number between 0
185
Cu is the harvest cost of the uth harvest, C′u the is cost of cleaning up and reforesting the
186
denotes the set of admissible thinning at stand age tu.
184
insect-attacked stand area, C0 is establishment and silvicultural costs, r is interest rate, σtu
187
188
Gross harvest revenues of the uth harvest are summed over n trees, j timber assortments as
189
products of prices pj (dollars per cubic metre), harvested volumes gj (m3 ha−1), and salvaged
190
volumes s ∙ Lj (m3 ha−1 ). The net harvest revenues are attained by subtracting harvest costs,
192
Cu and C′u , consisting of felling and on-site transport costs. The net harvest revenues are
193
interest.
191
discounted to the beginning of the rotation period by a factor (1 + r)−tk , where r is the rate of
194
195
The bare land value in Equation 2 is maximised subject to the stand dynamics that define the
196
stand state just prior to the harvest at age tu, and the catastrophe scenarios between present (tu)
8
197
and previous (tu-1) harvest. This is given as a function of stand state prior to the previous
198
harvest (Ztu−1 ), the intensity of the previous harvest (hu-1), the possibility of a catastrophe
199
200
happening and its damage intensity (btu−1 ) and the time difference between present tu and tu-1.
In our study, stand structure is described with n trees (n = 3), each tree is characterised by
201
variables reflecting its current dimensions (age, diameter, height etc.) and an expansion factor
202
representing the number of trees per hectare. All remaining trees are removed in the final (kth)
203
harvest.
204
205
The thinning rate equation 2.3 defines thinning rates for each diameter class. The thinning
206
rate is a piecewise linear function of tree diameter in our study, relative to the smallest and
207
largest stand diameters. The thinning parameters define the thinning rates at the corner points
208
of the piecewise linear function. The thinning definition can be adjusted to simulate different
209
thinning types by changing only a few parameters per thinning (for more details, see Valsta,
210
1992a).
211
212
2.4 Method for finding the optimal solution
213
214
Optimizing the management schedule means finding the optimal values for a set of decision
215
variables (DV), i.e. the optimal DV values for each thinning and final harvest (Miina, 1996;
216
Pukkala and Miina, 1997, 1998; Go´mez et.al, 2006). Because the number of thinnings is
217
not a continuous variable, schedules with different numbers of thinnings are to be treated as
218
separate optimisation problems. The management regime was specified by the number of
9
219
thinnings and the DVs, chosen as follows:
220
221
For each thinning: (a) years since previous thinning (which were set as a multiple of five
222
years from the previous thinning), (b) remaining basal area, and (c) harvest percentage.
223
224
For final cuttings: (a) years from the last thinning to the regenerative cut. Noticeably, the
225
criteria for grouping trees in our study are based on their DBH. Tree volumes and stumpage
226
values are thereby computed by diameter classes classified into small stems (with diameters
227
ranging from 4 cm to 13 cm), medium stems (with diameters ranging from 14 cm to 21 cm)
228
and large stems (with diameters great than 21 cm).
229
230
To determine the optimal management schedule, an initial solution (a supposition of the
231
optimal DV combination) was first fed into the computation. By doing this, the value of the
232
objective function with these DV values and the given set of parameters (timber prices, unit
233
costs and discounting rate) was calculated. The vector of DVs consists of times between
234
forest operations (adjusted by 1 year) and information on how the operations is executed, e.g.,
235
thinning percentages at different tree diameters (adjusted by 0.1%). Based on the feedback
236
from the simulation system, changes in the DV values were made in an effort to improve the
237
schedule, and the objective function was once again calculated using the new DV values. The
238
optimal management schedule for a given number of thinnings was eventually found, after
239
this search-process several times was repeated several times once a defined convergence
240
criterion was met. The optimisation was carried out for 0, 1, 2, 3, 4 and 5 thinnings separately,
10
241
in order to find the number of thinnings that maximised the bare land value. In order to gain a
242
better economic return, thinning intensity limits for small and medium trees were additionally
243
set at 40%, which a previous study found to be effective for larch plantations in northeastern
244
China (Sun et al., 2005).
245
246
3. Results
247
248
3.1 Deterministic optimal solutions at 3.5% interest rate
249
250
Table 1 shows the deterministic optimal solutions at a 3.5% interest rate for different plots.
251
The symbols in the tables are: DBH is diameter at breast height (cm), S-trees (%) denotes the
252
portion of small trees removed, M-trees (%) denotes the portion of medium trees removed,
253
L-trees (%) denotes the portion of large trees removed, Ro.dg denotes mean diameter at the
254
end of rotation (cm), M.A.I denotes mean annual increment in cubic metres per year, SI
255
denotes site index, which is the dominant height in meters at the index age.
256
257
Computations of sawlog and pulpwood amounts were based on the diameter classification.
258
For computing the amounts of pulpwood in our study, inferior trees such as small stems and
259
dead or damaged trees (natural mortality or caused by pest outbreak) were used as pulpwood;
260
the computation of sawlog amounts are based on the volumes of medium and large stems.
261
262
Table 1
11
263
264
Table 1 shows the proportion of trees removed in optimal thinnings for all plots at a 3.5%
265
interest rate. The thinning type is determined with three variables that define the thinning
266
rates for minimum, middle and maximum diameters. The thinning rates for other diameters
267
are interpolated.
268
269
To maximise the objective function according to various sample plots, two to five thinnings
270
were required. Two thinning were required on Yiershi 48, three thinnings were required on
271
Dulaer 76, Yiershi 53, Yiershi 57 O-2 and Yiershi 57 O-3, four thinnings were required on
272
Dulaer 80-2, Yiershi 58 and Yiershi 57 O-1 and five thinnings were required on Yiershi 57.
273
The sample plots in our study varied in their initial stand age, site fertility, initial stand
274
density and geographic location. Prior management treatments also possibly varied. However,
275
bare land values in our study were calculated assuming no prior harvest revenues or previous
276
establishment costs.
277
278
3.2 Sensitivity analysis of varying numbers of thinning
279
280
Deciding on the number of thinnings and the rotation length constitute seperate problems with
281
corresponding numbers of variables. The number of thinnings is therefore a parameter given
282
to the simulation-optimisation system at the outset. The maximum bare land value for
283
different numbers of thinnings (i.e. 0, 1, 2, 3, 4 and 5 thinnings) on various plots under
284
optimal rotation lengths are shown in Figure 2. Although bare land value increased with the
12
285
number of thinnings until the highest was reached, most of the economic gain was achieved
286
with just one thinning. Four thinnings were practically as good as the optimal number (i.e. 3
287
thinnings) for plots such as Dulaer 76 and Yiershi 53. Five thinnings were practically as good
288
as the optimal number (i.e., 4 thinnings) for Yiershi 58 and Yiershi 57 O-1. For other plots,
289
the bare land values change limitedly with the number of thinnings varying by one.
290
291
Figure 2.
292
293
3.3 Effects of catastrophic insect outbreak
294
295
To reveal the effects of catastrophic events, or pest outbreaks in our study, scenarios with
296
different probabilities and the intensity of occurring catastrophes were applied. We used 100
297
scenarios to gain a better picture of the result variability, as the application of fewer scenarios
298
may cause instability of the results.
299
300
In even aged management, one key decision variable is rotation age (Buongiorno and Gilless,
301
2003). The same method of scenario approach was used to produce good statistical statements
302
regarding the effect of the rotation age. For this purpose, 100 replications of the simulation at
303
different rotation ages were made, holding all parameters constant except for the string of
304
random numbers (e.g. catastrophe timing and intensity). These 100 scenarios of insect
305
damage were then generated randomly by using the Poisson process and used to simulate the
306
distribution of present value of each of the optimal regimes obtained from the deterministic
13
307
analysis. The largest and smallest bare land values, as well as the mean bare land value, the
308
standard deviation of the bare land value and the standard error of the mean, were observed in
309
these 100 replications with different scenarios for each different rotation. Thes results are
310
gathered and showed in Table 2. Table 2 shows that in these 100 replications, the bare land
311
value ranged from a minimum value to a maximum value. For example, at rotation age 49,
312
the bare land value for plot Dulaer 76 ranged from a minimum of $6749.0 ha-1 to a maximum
313
of $12850.8 ha-1. The mean bare land value was $10218.6 ha-1, with a standard error of
314
$131.2 ha-1. The 95% confidence interval of the mean bare land value is: 10218.6 ± 2 × 131.2
315
= ($9956.2 ha-1, $10481.0 ha-1). This 95% confidence interval does not contain the bare land
316
value obtained through the deterministic simulation, $11989.5 ha-1. Thus, the mean results
317
from the stochastic simulation are very different from those of the deterministic simulation.
318
319
Also, the catastrophes have an effect on the numbers of thinning when seeking a highest bare
320
land value. Table 3 indicates that the numbers of thinnings for most plots are decreased in the
321
stochastic case.
322
323
Table 3.
324
325
In Table 3, number of thinnings was decreased from four to three for Dulaer 80-2, Yiershi 58
326
and Yiershi 57 O-1; for Yiershi 53 and Yiershi 57 O-3, number of thinnings was decreased
327
from three to two; for Yiershi 57, number of thinnings was decreased from five to four.
328
14
329
Details of Table 2 are displayed in the Appendix. The optimal solution for each plot in the
330
deterministic analysis was added to perform the comparisons. The optimal bare land value for
331
both deterministic and stochastic cases are in bold characters.
332
333
Table 2 indicates that the mean bare land values from the stochastic simulation (with 100
334
replications of different scenarios) are 14.8% to 22.9% lower than the deterministic
335
simulation, which corresponds to the different plots.
336
337
When the stand was subject to the risk of catastrophe, increasing risk reduced the optimum
338
rotation. This is demonstrated in Figure 3, we took Dulaer 76 as example, where the
339
deterministic case is compared with three levels of risk, presented by annual probabilities of a
340
catastrophic insect outbreak. Damaged trees are logged by thinnings.
341
342
Figure 3.
343
344
Given a constant thinning rate, the optimum rotation was shortened from 53 yr to 47 yr when
345
the probability increased from 0 to 0.04. Thinning was brought forward in higher risk cases
346
(Table 4).
347
348
Table 4.
349
350
The other decision criterion examined in this study was conditioned by a specified bare land
15
351
value. By adjusting rotation length, the forest manager wishes to maximize the probability of
352
attaining a target bare land value. In this formulation, the decision-maker can increase
353
willingness to take risks by raising the bare land value requirement. When thinning rate was
354
constant, increased risk-taking caused the optimum rotation, i.e., the one with highest
355
probability, to decrease (Figure 4). For example, the probability that a simulated bare land
356
value is high than 10000 USD was 0.63 when the rotation was 49 years. Namely, in 63 of the
357
100 simulations made, the bare land value was higher than 10000 USD.
358
359
Figure 4.
360
361
4. Discussion
362
363
The main purpose of the present study is to reveal the effects of insect outbreaks on the
364
financial performance of Dahurian larch. To do this, a stand level optimisation approach is
365
applied. Stochasticity is introduced in the form of scenarios.
366
367
Results in the deterministic case indicate that stands with high mean annual increment such as
368
Yiershi 48, Yiershi 57 O-2 and Yiershi 57 O-3, produced the highest bare land values among
369
all the plots. Many factors affect rotation age, e.g. planting density and thinning frequency
370
and intensity (Hyytiäinen, 2005; Go´mez et.al, 2006).
371
372
Clearcutting would take place around the age of 42 to 58 years on various plots. The mean
16
373
diameters of all the plots at the end of the rotation period in this deterministic case ranged
374
from 15.9 cm to 21.3 cm, which is different to results received in other studies (e.g. Li, 1987;
375
Sun et al., 2001; Shi and Feng, 2005). These studies had longer rotation ages and greater
376
mean diameters (usually > 20 cm) at the end of rotation period. On the other hand, Sun et al.
377
(2005) and Gao et al. (2009) reported similar results of mean diameters at the end of rotation
378
in their research as are reported by our study. The heavy cutting intensity of large trees in
379
previous thinnings and the variance of the initial state and site quality might be key reasons
380
for the difference in mean diameters at the end of rotation. The optimal rotation period for
381
plot Yiershi 57 O-3 is much shorter compared to other plots. Usually larch plantation stands
382
must be 45 years or older to meet clearcutting requirements (Gao et al., 2009). However,
383
some earlier studies reported that the optimal rotation can be shorter than recommended. For
384
example, Sun et. al (2005) claimed that the rotation period of larch is 37–39 years with the
385
first thinning timed 25 years. This result was naturally found using a different discount rate,
386
timber price and site location.
387
388
For the stochastic analysis, the risk of catastrophic pest outbreaks shortens the optimum
389
rotation considerably and decreases the numbers of thinning, respectively. Among all nine
390
plots, rotation period of Yiershi 53 decreased mostly in stochastic case, by eight years. For
391
Yiershi 57 the risk of catastrophic pest outbreak shortened the optimum rotation period by
392
seven years, for Dulaer 80-2 and Yiershi 57 O-1 both by six years. And for most plots,
393
numbers of thinnings were reduced by once than the deterministic optima. The increasing
394
catastrophic risks affected the bare land values with variability. The bare land values from the
17
395
stochastic case were approximately 14.8% to 22.9% lower on average than in the
396
deterministic case. Compare to the deterministic case, bare land value of Yiershi 57 had the
397
largest decreasing of 22.9% lower, and bare land value of Dulaer 76 had the smallest
398
decreasing of 14.8% lower among all nine plots. In addition, with simulating a specified
399
probability (%) we could achieve at least as good a bare land value as with the deterministic
400
case when e.g. Dulaer 49-yr rotation period is applied in the stochastic case (deterministic
401
case 53-yr rotation).
402
403
Our results point out that initial conditions and management plans may strongly influence the
404
results of including stochasticity into stand level optimisation. The source of stochasticity was
405
the occurrence of a catastrophic pest outbreak. To determine the solution with confidence and
406
get a better picture of the variability of the results, 100 scenarios were exerted—a smaller
407
number of scenarios application may cause instability and report incorrect solution.
408
409
Information about possible outcome variability has additionally been gained from the
410
stochastic simulation. While the law of averages may be relevant for owners with several
411
stands of this type, such as the owner in our study, it may not be relevant for owners with a
412
single woodlot. Any one of the stochastic simulation outcomes is possible for this owner
413
group, and outcome variability will affect property values and influence managerial decisions.
414
415
In our study, the scenarios were formed by generating random outcomes of the defined
416
stochastic process (i.e. catastrophic pest outbreak). Another alternative would be to form
18
417
scenarios subjectively, so that they would be likely realizations of the future. Furthermore,
418
one of the strengths of the present scenario-based stochastic simulation is that the distribution
419
of the objective function values is known and can be analysed with respect to relevant
420
variables, by scenario approach. The scenario simulations produce huge amounts of
421
information which can be used to evaluate management strategies.
422
423
Quantitative probabilistic risk analyses in our study used mathematical definitions that
424
described the relationship between the expected value of the probability of the event
425
occurring at a particular location and intensity, and the net value change given that the event
426
has occurred. Calculating the probabilities for stochastic processes and the effects of
427
catastrophic or expected loss in this manner is advantageous for risk analysis related to
428
disturbance processes. The results and approach offer managers several opportunities because
429
they provide information that can be used to assess risks that are difficult to amend through
430
mitigation or in situations where additional mitigation investments have limited effects.
431
432
Forest planning and decision-making usually deal with large forest area, instead of a single
433
stand (Haight and Monserud, 1990; Valsta, 1992b). The applicability of the presented system
434
depends on the production objectives of the forest. If the objective is to maximise net income
435
or profitability while considering a stochastic situation or risk of catastrophe, the system can
436
be used directly. It is possible to separately optimise the management of each stand, or use the
437
system to produce management instructions for different sites and stand densities in similar
438
circumstances. Stand-level optimisations also serve comparative analyses on the effects of
19
439
economic or biological factors on stand management (Valsta 1992b).
440
441
We believe that the system presented in our study is a useful contribution towards improved
442
decision-making in the management of Dahurian larch stands under risk of catastrophic insect
443
outbreaks. The system enables us to assess the likelihood and severity of economic
444
consequences when exposed to an invasive species. The system is essentially a stand level
445
decision-support tool, which can be used to simulate alternative management regimes and
446
find economic management strategies. The system can also be used in forest-level
447
optimisation in conjunction with linear programming or other optimisation procedures.
448
System reliability will obviously depend on the reliability of the simulation system. It is
449
therefore important to refine and enhance it, e.g. by more precisely quantifying economic
450
and ecological damage related to the invasive species; incorporating stochastic structure into
451
the simulation system (Fox et al. 2001) which leads to valid statistical inference, improved
452
estimation efficiency and more realistic and theoretically sound prediction. It is proposed that
453
individual-tree modelling methodologies need to characterise and include structured
454
stochasticity in future research.
455
456
Funding
457
458
This work was supported by the Academy of Finland [grant number 114201, Biodiversity and
459
Forest Pest Problems in Northeast China, BIOPROC].
460
20
461
Acknowledgements
462
463
We wish to thank the helpful staff of the Aershan Forestry Bureau in Inner Mongolia,
464
Aershan, P.R China.
465
466
References
467
468
469
Buongiorno, J. and Gilless, J. K. 2003. Decision methods for forest resource management.
CA: Academic Press, San Diego, pp. 313–336.
470
Fang, Z., Borders, B.E. and Bailey, R.L. 2000. Compatible Volume-Taper Models for
471
Loblolly and Slash Pine Based on a System with Segmented-Stem Form Factors. Forest Sci.
472
46(1), 1–12.
473
474
FAO-UNESCO. 1989. Soil Map of The World (revised legend). Tech. Pap. 20, ISRIC,
Wageningen, Netherlands, 138 p.
475
Forest Diseases and Insect Pests Control and Quarantine Station. 2002. Internal
476
Communication. Aershan Forestry Bureau of Inner Mongolia, Aershan 137801, P.R China.
477
Forest Diseases and Insect Pests Control and Quarantine Station. 2006. Meteorological and
478
Hydrological Records (Internal Communication). Aershan Forestry Bureau of Inner
479
Mongolia, Aershan 137801, P.R China.
480
Forsell, N. 2003. Planning Under Risk and Uncertainty: Optimizing Spatial Forest
481
Management Strategies. Doctoral Thesis, Swedish University of Agricultural Sciences,
482
Umeå, Sweden.
21
483
484
Fox, J.C., Ades, P.K. and Bi, H. 2001. Stochastic structure and individual-tree growth models.
Forest Ecol. Manage. 154, 261–276.
485
Gao, H., Li, F. and Jia, W. 2009. Study on Diameter Class Distribution and Economic
486
Rotation of Larix gmelini Plantation. Forest Eng. Vol.25, No.4, 10–14. (In Chinese, with
487
English summary)
488
T. Go´mez , M. Herna´ndez , M.A. Leo´n, R. Caballero. 2006. A forest planning problem
489
solved via a linear fractional goal programming model. Forest Ecol. Manage. 227, 79-88.
490
González-Olabarria, JR. and Pukkala, T. 2011. Integrating fire risk considerations in
491
492
493
landscape-level forest planning. Forest Ecol. Manage. 261, 278–287.
Haight, R.G. and Monserud, R.A., 1990. Optimising any-aged management of mixed-species
stands. II. Effects of decision criteria. Forest Sci. 36 (1), 125–144.
494
Hanewinkel, M., Peltola, H., Soares, P. and González-Olabarria, J.R. 2010. Recent
495
approaches to model the risk of storm and fire to European forests and their integration into
496
simulation and decision support tools. Forest Syst. 19(SI), 30 –47.
497
498
499
500
Huber, M. 2010. Assessing the long-term species composition predicted by PROGNAUS.
Forest Ecol. Manage. 259, 614-623.
Hyytiäinen, K., Tahvonen, O. and Valsta, L. 2005. Optimum Juvenile Density, Harvesting,
and Stand Structure in Even-Aged Scots Pine Stands. Forest Sci. 51(2), 120 –133.
501
Jenkins, M.J., Hebertson, E., Page, W. and Jorgensen, C.A. 2008. Bark beetles, fuels, fires
502
and implications for forest management in the Intermountain West. Forest Ecol. Manage.
503
254, 16–34.
504
Jiang, S. 2009. Study on Forest Dynamic Model for Larix in Plantation. Doctoral Thesis,
22
505
Northeast Forestry University, Hei Longjiang, China. (In Chinese, with English summary)
506
Li, F. 1987. Study on Diameter Distribution and Models Predicting Yields for Natural
507
Dahurian Larch stands. J. NE Forestry Uni., Vol.15, No.4, 8–16. (In Chinese, with English
508
summary)
509
510
511
512
Miina, J. 1996. Optimizing thinning and rotation in a stand of Pinus sylvestris on a drained
peatland site. Scan. J. For. Res. 11, 182–192.
Mitchell, S.J. 2013. Wind as a natural disturbance agent in forests: a synthesis. Forestry 86
(2), 147–157.
513
Monserud, R.A. and Sterba H. 1996. A basal area increment model for individual tree
514
growing in even- and uneven-aged forest stands in Austria. Forest Ecol. Manage. 80,
515
57–80.
516
517
518
519
Monserud, R.A. and Sterba H. 1999. Modeling individual tree mortality for Austrian forest
species. Forest Ecol. Manage. 113, 109–123.
Moore, H.T., Van Miegroet, H. and Nicholas, N.S. 2008. Examination of forest recovery
scenarios in a southern Appalachian Picea–Abies forest. Forestry 81, 183–194.
520
Peltola, H., Kellomäki, S., Väisänen, H. and Ikonen, V.P. (1999). A mechanistic model for
521
assessing the risk of wind and snow damage to single trees and stands of Scots pine,
522
Norway spruce, and birch. Can. J. Forest Res. 29(6), 647–661.
523
524
525
526
Pukkala, T. and Miina, J. 1997. A method for stochastic multiobjective optimisation of stand
management. Forest Ecol. Manage. 98, 189–203.
Pukkala, T. and Miina, J. 1998. Tree-selection algorithms for optimizing thinning using a
distance-dependent growth model. Can. J. For. Res. 28, 693–702.
23
527
528
529
530
Reed, W.J. 1984. The effects of risk of fire on the optimal rotation of a forest. Journal of
Environmental Economics and Management 11:180-190.
Reed, W.J. and Errico, D. 1985. Assessing the long-run yield of a forest stand subject to the
risk of fire. Can. J. For. Res. 15, 680-697.
531
Shi, L.P. and Feng, Z.K. 2005. Basic methods for establishing plantation stand growth
532
estimate models. J. Beijing Forestry Uni., Vol. 27, Supp. 2, 222–225. (In Chinese, with
533
English summary)
534
Sun, C., Shen, G., Li, J. and Jia, L. 2001. Study on the Present Condition and the
535
Potentialities of the Productivity of Main Tree Species Plantation of China. Forest Res.
536
14(6), 657–667. (In Chinese, with English summary)
537
Sun, J., Li, S. and Wang, Y. 2005. Consideration of Short Regulatory Rotation Age of Larch
538
Plantation. Forestry Sci. Technol. Inf., Vol. 37, No. 2, 4–5. (In Chinese, with English
539
summary)
540
541
542
543
544
545
Valsta, L. 1992a. A scenario approach to stochastic anticipatory optimization in stand
management. Forest Sci. 38, 430–447.
Valsta, L., 1992b. An optimisation model for Norway spruce management based on
individual-tree growth models. Acta For. Fennica 232, 20p.
Wang, L., Luo Y., Huang, H., Shi, J., Heliövaara, K., Teng, W. et al. 2009. Reflectance
Features of Water Stressed Larix gmelinii Needles. For. Stud. China 11(1), 28–33.
546
Wang, X., He, H.S., Li, X., Chang, Y., Hu, Y., Xu, C. et al. 2006. Simulating the effects of
547
reforestation on a large catastrophic fire burned landscape in Northeastern China. Forest
548
Ecol. Manage. 225, 82–93.
24
549
550
Appendix
551
552
The appendix gives a more detailed description of the simulator constructed for our current
553
study, based on Monserud and Sterba (1996) and Jiang (2009). The growth models were
554
estimated for periods of five years. When shorter growth periods are needed, they are linearly
555
interpolated from the full 5-year periods.
556
557
The simulation starts with a tree list of Dahurian larch containing data on tree DBH, height,
558
biological age, age at breast height, and the number of trees per hectare. Before the growth
559
projection, a set of stand variables is computed including stand basal area, an individual tree’s
560
percentile in the stand basal area distribution, and mean height. The simulation of one 5-year
561
time span consists of the following four steps:
562
563
1.
Tree basal area increment:
564
565
ln(BAI) = 0.7144 + b SIZE + c COMP + s SITE,
(E.1)
566
567
where BAI is the 5 year basal area increment, SIZE stands for the tree size variables, COMP
568
stands for the competition variables, and SITE stands for the site variables.
569
570
In detail, diameter at breast height was used as a measure of stem size and crown ratio as a
25
571
measure of crown size, yielding the expression:
572
573
b SIZE = 1.1547 ln(DBH) – 0.000101 DBH2 + 0.4999 In(CR),
(E.1.1)
574
575
where DBH (cm) is the diameter at breast height and CR is the crown ratio. Although DBH
576
and CR here are referred to as size variables, they also reflect the effects of past competition.
577
578
Current competition effects were expressed as:
579
580
c COMP = – 0.0159 BAL,
(E.1.2)
581
582
where BAL is the basal area (m2 ha-1) of trees larger in diameter than the subject tree. The
583
BAL of the largest tree is 0.0, and for the smallest tree equals the stand basal area minus that
584
tree's basal area.
585
586
Site description is very important as it influences tree growth, site function was expressed as:
587
588
s SITE = – 0.00174 (ELEV)2 + 0.1882 S + 0.3309 V,
(E.1.3)
589
590
where ELEV is the elevation in hectometers (100 m), S is a set of dummy variables for
591
various soil groups (the soil group of our study site is Spodi-Dystric Cambisol on silicate
592
material); V is a set of dummy variables for various vegetation types (the vegetation type of
26
593
our study site is Oxalis acetosella types (Moder humus in conifer stands)). Here, the values
594
for the corresponding dummy variables are S = 1 and V = 1. (for more details, see
595
FAO-UNESCO, 1989; Jiang, 2009; Forest Diseases and Insect Pests Control and Quarantine
596
Station, 2002)
597
598
2. Live crown ratio:
599
600
ln(CR) = – 0.95042 – 0.008886239175 BA + 0.30066 ln(DBH) – 0.59302 ln(H) + 0.19558
601
ln(PCT),
(E.2)
602
603
where BA is stand basal area (m2 ha-1), DBH is diameter at breast height (cm), H is tree
604
height (m), and PCT is a tree’s percentile in the stand basal area distribution.
605
606
3.
Height growth curve:
607
609
H= 1.3 + 1.3432 Ab 0.7571 [1 – exp (0.02296 SI 0.7548 DBH)],
610
where H is tree height (m), Ab is age at breast height, SI is site index, DBH is diameter at
611
breast height (cm). Additionally, the site index in Equation (E.3) was decided using the
612
following equation:
608
(E.3)
613
614
SI = Hmean * exp(
12.51004
A
–
12.51004
A1
),
(E.3.1)
27
615
616
where Hmean is the stand mean height (m), A is the stand age, A1 is the index age (the index
617
age for larch plantations is 30 years).
618
619
4.
Tree mortality:
620
621
P = [1 + exp (4.407 – 12.9395 / DBH + 2.2039 CR – 0.0326 BAL)] -1，
(E.4)
622
623
where P is the probability of mortality (in a 5-year period), DBH is diameter at breast height
624
(cm), CR is crown ratio, BAL is the basal area of larger trees (m2 ha-1), and b0 - b5 are
625
parameters. Furthermore, DBH, CR and BAL are values at the beginning of the 5-year period.
626
627
The merchantable volumes of individual trees are predicted and computed using the
628
segmented stem taper and volume equations modified from Fang et al. (2000):
629
630
631
632
633
d2= c2 H
k−b1
b1
(1 − z)
k−β
β
α1 I1 + I2 α2 I2 ,
(E.5.1)
k
k
Vm = c2 H b1 [b1 + (I1 + I2) (b2 – b1) t1 + I2 (b3 – b2) α1 t 2 – β (1 − z) β α1 I1 + I2 α2 I2 ],
(E.5.2)
634
635
636
where
2
c
=
I1 = �
1, 0.076907 ≤ z ≤ 0.56896
,
0,
otherwise
d3 −
k
b1
d1 DBHd2 H
b1 (1− t1 ) + b2 (t1 − α1 t2 ) + b3 α1 t2
I2 = �
1, 0.56896 < z ≤ 1
,
0,
otherwise
k
k
, t1 = 0.923093b1 , t 2 = 0.43104b2 , α1 =
28
637
638
0.923093
(b2 −b1 )k
b1 b2
, α2 = 0.43104
(b3 −b2 )k
b2 b3
, β= b11−(I1+I2) b2 I1 b3 I2 , z =
h
H
, k =
π
40000
,
b1 = 0.000010481, b2 = 0.000038139, b3 = 0.000027594, d1 = 0.000019913, d2 = 1.94675, d3
640
= 1.185269, Vm is the merchantable volume (m3), d (cm) is the diameter of the measurement
641
aboveground height to the measurement point on the stem.
639
point at height h (m), DBH (cm) is diameter at breast height, H is tree height (m), h is the
642
643
644
645
646
z = 0.076907 and 0.56896 are the relative locations of the two inflection points of Dahurian
larch trees, which divided a stem into three corresponding stem-segments. A special case
included in the merchantable volume equation occurs, when h = H, Vm arises as the total
stem volume.
647
648
Log volumes obtained from the model are based on stem dimensions. The minimum sawlog
649
diameter and length are 14 cm and 2m, respectively; the minimum pulpwood diameter and
650
length are 5 cm and 2m, respectively.
651
h
653
From Equation E.5.1, it is possible to calculate the value of z (z =
654
measurement point where d is a specific value (e.g. d = 14 cm). Bringing this value z into
Equation E.5.2, we can compute the sawlog volume (with diameter ≥ 14 cm). Using the same
655
method, we can also calculate the pulpwood volume (with diameter ≥ 5 cm and < 14 cm).
652
656
657
Table 2.
29
H
) at a specific