Journal of Forest ScienceVol. 24, No. 3, pp. 143~149, December 2008

Economic Analysis of Snow Damage on Sugi
(Cryptomeria japonica) Forest Stands in Japan Within theForest Stand Optimization Framework
Atsushi Yoshimoto1, Akio Kato2, and Hirokazu Yanagihara3*
1

Dept. of Math. Analy. and Stats. Infer., Institute of Statistical Mathematics, Tokyo, Japan, yoshimoa@ism.ac.jp
Toyama Prefectural Agricultural, Forestry & Fisheries Research Center, Toyama, Japan, kato-a@fes.pref.toyama.jp
3
Dept. of Mathematics, Hiroshima University, Hiroshima, Japan, yanagi@math.sci.hiroshima-u.ac.jp

2

ABSTRACT : We conduct economic analysis of the snow damage on sugi (Cryptomeria japonica) forest stands in Toyama Prefecture, Japan. We utilize a single tree and distant independent growth simulator called “Silv-Forest.” With this growth simulator,
we developed an optimization model by dynamic programming, called DP-Silv (Dynamic Programming Silv-Forest). The MSPATH (multiple stage projection alternative technique) algorithm was embedded as a searching algorithm of dynamic programming.
The height / DBH ratio was used to constrain the thinning regime for snow damage protection. The optimal rotation age turned
out to be 65 years for the non-restricted case, while it was 50 years for the restricted case. The difference in NPV of these twocases as the induced costs ranged from 179,867 to 1,910,713yen/ha over the rotation age of 20 to 75 years. Under the optimalrotation of 65 years, the cost became 914,226 yen/ha. The estimated annual payment based on the difference in NPV, was from9,869 yen/ha/yr to 85,900 yen/ha/yr. All in all, 10,000 yen/ha/yr to 20,000 yen/ha/yr seems to cover the payment from the rotation
age of 35 to 75 years.
Keywords : Snow damage, Cryptomeria japonica, Algorithm, Forest management

INTRODUCTION

snow damage on sugi forest stands in Toyama Prefecture,
Japan. We utilize a single tree and distant independent

Snow damage on artificial forest stands has been one

growth simulator called “Silv-Forest.” This growth simulatorof the main concerns in snow regions, Japan. It breaks

was originally developed by Tanaka (1995) and modified

stems and branches with the large economic loss to forest

by Kato and Tanaka (1997) and Kato et al. (2008) for

owners. Kato et al. (1992) conducted the analysis on thesugi forests in Toyama. One of the main characteristics ofratio of height to DBH (diameter at breast height) (hereafter

Silv-Forest is that the DBH distribution is first createdcalled H/D ratio) for a sugi (Cryptomeria japonica) forestfrom the observed data, then the distribution progresses

stand against snow damage, and found out that the H/Dmainly based on diameter growth at each DBH classratio should be maintained less than 65% to avoid thefollowing the normal distribution with the derived mean

damage under the thinning regime. Although this require-of growth and the given coefficient of variation for growth.

ment seems to reduce the benefits from the management,

It has been used on the spreadsheet basis of Microsoft

it could be offset by the benefits of avoiding the damage

Excel Software. With some arrangements, we modify it

beforehand.

further to be optimized for the thinning regime after

In this research, we conduct economic analysis of theconverting it into the FORTRAN programming language.

* Corresponding author: (E-mail) yoshimoa@ism.ac.jp※ Advance in Forest Management and Inventory-selected papers from the international conference of IUFRO (Div. 4.01, 4.02, 4.04),
Chuncheon, Korea, Oct. 13-17, 2008.

144 ‧ Journal of Forest Science

With this growth simulator, we develop an optimization

MS-PATH (Multiple Stage PATH) algorithm to overcome

model by dynamic programming, called DP-Silv (Dynamic

the problem by considering all possible future path com-

Programming Silv-Forest). The MS-PATH (multiple stage

bination. The recent application of MS-PATH can beprojection alternative technique) algorithm is embedded as

found in Yoshimoto and Marušák (2007). In what follows,a searching algorithm of dynamic programming. The H/D

we elaborate the algorithm.

ratio is used to constrain the thinning regime for snow

Let x(t) be a state variable representing a forest standdamage protection. The paper is organized as follows. In

at time t, and ẋ (t) be its first derivative with respect tothe second section, the MS-PATH algorithm is mathemati-time t. Introducing a marginal function, π(․), of the objec-

cally elaborated. In the third section, we conduct econo-

tive for the management, say, the net present value of the

mic analysis on the snow damage, then some concluding

profits, over time, the objective is to maximize its integral

remarks are provided in the last section.

from the initial state (plantation), x0 at time t0, to the ending state (final harvest) xn at time tn with respect to a

Dynamic Programming Specification
The dynamic programming approach has been developedand extensively applied for optimizing thinning regimes

control variable, ẋ (t).
tn

max J = ∫ π ( x(t ), x& (t ), t ) dt
{ x& ( t )}t0

x(t0 ) = x0

(Arimizu, 1958a, b; Amidon and Akin, 1968; Schreuder,(1)

x(tn ) = xn

1971; Kilkki and Väisänen, 1970; Adams and Ek, 1974;
Brodie et al., 1978; Brodie and Kao, 1979; Chen et al.,In the thinning management problem, thinning inten-

1980; Martin and Ek, 1981; Riitters et al., 1982; Brodiesity, T(t) at time t is of the major concern, and controlsand Haight, 1985; Haight et al., 1985; Valsta and Brodie,

x(t) and ẋ (t), so that the control variable can be replaced

1985; Arthand and Klemperer, 1988; Torres and Brodie,

by thinning, T(t),1990). However, as forest stand growth models weredeveloped with growth complexity, the dynamic programming approach became inefficient for optimization due totn

max J = ∫ π ( x(t ), T (t )) dt{T ( t )}t0

x(t0 ) = x0

vast increase in dimensionality for the state description of(2)x(tn ) = xn

iteration called the curse of dimensionality.
While the curse of dimensionality had been the main

Converting the objective function into the discrete

problem for these early applications of dynamic pro-thinning problem with thinning at age of t (t = t0, t1, …,

gramming. Paredes and Brodie (1987) introduced a newtn), we have,searching algorithm called PATH (Projection Alternative
Technique) within the framework of both network theory
and the theory of the Lagrange multipliers to resolve theproblem. Their derivation of the algorithm was somewhat
ad-hoc with difficulty of searching for optimal values of
Lagrange multipliers. With the concept of calculus ofvariation, Yoshimoto et al. (1988) introduced a new derivation of the algorithm, and pointed out the shortcomingof PATH in considering a long-term effect of thinning
activities on the objective function. They developed thenmax J = ∑ ∫ π ( x(t ) | T (ti −1 )) dt{T ( ti )}

i =1

ti

ti−1

x(t0 ) = x0

(3)x(tn ) = xnIntroducing the integrand of Π(x(ti)󰠐T(ti-1)) =  π(x(τ)󰠐


T(ti-1))dτ, the objective function of equation becomes,

Economic Analysis of Snow Damage on Sugi (Cryptomeria japonica) Forest Stands in Japan Within the … ‧ 145nstage for two sequential thinning over multiple stages.

max J = ∑ ∫ π ( x(t ) | T (ti −1 )) dt{T ( ti )}

i =1

ti

Note that Ti,i-j is thinning intensity at time ti-j with the

ti −1

n

= ∑ {Π ( x(ti ) | T (ti −1 )) − ∫ π ( x(t ) | T (ti −1 )) dt} (4)
ti−1

0

i =1

T

elapse of the stage j from the i-th stage, Πi(Ti,i-j) is thecontribution value of a forest stand at time ti after having

Let us introduce the contribution value function, Π (T(ti)),thinning Ti,i-j at time ti-j, ΠT(Ti,i-j) is the contribution valueof thinning, T(ti), which satisfies the following,of thinning Ti,i-j implemented at time ti-j, j* is an optimal*ti∫ π ( x(t ) | T (t )) dt = Π ( x(t ) | T (t )) − Π (T (t ))
T

i

0

i −1

i

i

elapse of the stage from the i-th stage, T i,i-j is an optimal(5)

The first term of the right-hand side is the contribution
value of a forest stand, Π(․), at time ti after having thinning T(ti-1) at time ti-1 and the second term is the amounttaken by thinning T(ti) at time ti. As a result, the difference is the contribution value of a stand just after having thinning at time ti. With equation (5), equation (4)becomes.

(T*i,i-j*), is the contribution value of a forest stand at time*

ti after having an optimal thinning intensity T i,i-j* at time
ti-j* with an optimal elapse of the stage j*. The algorithm
searches for an optimal solution by maximizing the sumof the contribution value of a forest stand and thinning.

Analysis
Regarding the growth simulator of Silv-Forest, we im-nmax J = ∑ {Π ( x(ti ) | T (ti −1 )) − ∫ π ( x(t ) | T (ti −1 )) dt}{T ( ti )}thinning intensity at time ti-j targeting time ti, Π*i = Π

ti−1

plemented the following modification. First we utilized the

0

i =1
n

= ∑ {Π ( x(ti ) | T (ti −1 )) + Π T (T (ti −1 )) − Π ( x(ti −1 ) | T (ti − 2 ))}
i =1

(6)

Note that only the first two terms on the right-handside have affects from thinning, T(ti).

relative taper curve function to estimate volume of a tree
by every 1/10 segmentation of a tree. Second, since theoriginal Silv-Forest did not consider self-thinning, we embedded it to the model.
The taper curve function is to estimate the i-10 th rela-

With one stage look-ahead search in the PATH

tive diameter size, Di to DBH, by the following empiricalalgorithm (Paredes & Brodie, 1987), the following

polynomial function of a relative height, xi, with all para-

optimality equation is obtained.

meters estimated as a function of the corresponding DBH

fi *= max{ f i (Ti −1 )}fi (Ti −1 )= Π i (Ti −1 ) + Π (Ti −1 ) − Π

and height, so that at each DBH class one unique taper
curve is derived.

{Ti−1 }

T

*
i −1

+f

*
i −1

(7)

where Ti = T(ti), Πi (Ti-1) = Π(x(ti)󰠐T(ti-1)), ΠT (Ti) = ΠT(T(ti)),
Πi* = Πi(T*i-1), and Ti* is an optimal thinning intensity attime ti. Considering influence of thinning over multiple
periods by MS-PATH, the optimality equation of MSPATH becomes,fi *i = 1,...,10(9)Each parameter is a function of DBH and height definedby the following bilinear form.
ai = b0i + b1i ⋅ H j + b2i ⋅ DBH j + b3i ⋅ H 2j + b4i ⋅ H j ⋅ DBH j += max { f i ,i − j (Ti ,i − j )}+ b5i ⋅ DBH 2j{Ti ,i − j , j }fi −i j (Ti ,i − j ) = Π i (Ti ,i − j ) + Π T (Ti ,i − j ) − Π *i − j + f i *− j

Di = a1 ⋅ xi + a2 ⋅ xi2 + a3 ⋅ xi3 + a4 ⋅ xi20

(8)

(10)

Converting these equations into the framework of theThis is to search for an optimal thinning intensity as

estimated generalized least square (EGLS) estimation (seewell as an optimal elapse of the stage, j, from the (i-j)-thVonesh and Carter, 1987), we have the following statisti-

146 ‧ Journal of Forest Sciencecal formulation for the j-th taper data.

above {b}a j : k × 1 vector of between-individual explanatory vari-

y j = X jβ j + ε j( j = 1,..., n),β j = Ξ′a j + ν jables of DBH and height of the j-th taper data,(11)νj : q × 1 between-individual error vector νj ~ i.d.E[νj] = 0,Cov[νj] = Δ.

This results in the following for the final estimation.
We used the data from the stem analysis to estimate
y j = X j Ξ′a j + X j ν j + ε j ( j = 1,..., n),coefficients. Figure 1 shows the fitted results of the(12)estimation with the results in Table 1.
The second modification is to incorporate a function of

where
yj : pj × 1 vector of repeated measurements from stem

self-thinning into the growth model. The following is theanalysis for diameter of the j-th taper data,

relationship of the number of trees and individual average

Xj : pj × q within-individual design matrix of relative

volume estimated with the use of the above relative taper
curve.

height of the j-th taper data,βj : q × 1 random regression coefficient of the j-thtaper data corresponding to the above {a},

1
1
v
=
−
N N 0 3470592 ⋅ N 0-0.9184

εj : pj × 1 within-individual error vector εj ~ i.d.E[εj],(13)Cov[εj] = σ2Ipj,
Ξ′ : k × q mean parameter matrix corresponding to thewhere N is the current number of trees, N0 is the initial
number of trees of plantation and v is the average individual tree volume at the current.
One variant of sugi called Kawaidani sugi was considered for the analysis. This is a typical variant in Toyama,
Japan. The stand age was 15 years old with average DBH
of 14.41 cm, average height of 9.39 m, average tree volumeof 0.079 m3 and TPH (trees/ha) of 2001. The main
growth component of Silv-Forest is that it deals with the
DBH distribution and its dynamics by the transition probability matrix derived from the induced normal distribution of a tree DBH growth. Figure 2 shows the DBH
distribution and its dynamics without any thinning over
55 year time horizon. The control variables for optimization were thinning timing, thinning intensity, and thinn-

Fig. 1. Fitted results of the relative taper curve.

Table 1. Estimated results for the relative taper curve












a1

0.937

0.516

-0.5920.004-0.032

0.034

a2

0.891

-0.9441.064-0.0590.190-0.144

a3

0.277

0.639

-0.9160.034-0.1170.099

a4-0.145-0.0420.117-0.0260.067-0.044

Economic Analysis of Snow Damage on Sugi (Cryptomeria japonica) Forest Stands in Japan Within the … ‧ 147

Fig. 3. NPV of the optimal solution.

Fig. 2. DBH distribution and dynamics.

diameter size, but also by the grade of log determined on
the basis of the position of log cut. The H/D ratio was

Table 2. Log price as a function of DBH for 4m size log

calculated by 100(H / D) by adjusting the dimension ofRange of Top Diameter
Grade of Log
of Log (cm)each estimate, where H and D are the average height (m)Price (Yen/m3)

and DBH (cm) of a forest stand.

10-13

0

9,714

We considered two cases for optimization. One was

14-18

1

9,792

without any H/D ratio requirement, and the other is with

14-18

2

9,792

20-28

1

18,94320-28

213,26430-38

123,46430-38

223,46440-1001

25,65940-1002

25,659Note: Grade of log is determined by the position of log
produced from a tree. “1” means the 1st log from thebottom of a tree, and the rest is the 2nd grade withsmall diameter as 0 grade

the ratio less than or equal to 65% at the thinning andfelling periods by

100 ⋅

Hn

H n ( m)≤ 65
DBH n (cm)(14)

: Average height at the n-th stage

DBHn : Average DBH at the n-th stage
The difference in the objective value from these twocases was regarded as a cost induced by the H/D ratio

ing method (thinning from below, thinning from above, or

requirement for snow damage protection. Figure 3 shows

proportional thinning).

the optimal net present value over different rotation age
3With 1% discount rate and 8,000 yen/m for fellingunder two cases. The gray “x” symbol represents a solu-

costs, we searched for an optimal thinning regime to maxi-

tion without the H/D ratio requirement, while the black

mize the net present value of all returns from thinning“x” is with the requirement. From these outcomes, theand final harvesting. Since our optimization model is

optimal rotation age turns out to be 65 years for the first

discrete, we applied 5 trees as an interval of thinning in-case, while 50 year rotation age was for the restrictedtensity and 5 years as an interval of thinning timing in thecase. The difference in NPV, CSD, ranged from 179,867

searching process of DP. The benefits were based on theto 1,910,713 yen/ha for the rotation age of 20 to 75 years.

price of logs produced from trees. We used the price data

Sudden change was observed for 25 and 30 year rotation.

on log as a function of the top diameter of log at the log

Under the optimal rotation of 65 years, the cost became

auction market in Toyama shown in Table 2. As can be914,226 yen/ha as in Figure 4.

seen from the table, the price changes not only by the top

The above difference was calculated based on different

148 ‧ Journal of Forest Science

20 years to 85,900 yen/ha/yr of the age 25 years. At the
optimal rotation age for the restricted case, it was 10,416
at age 50 years. More or less, 10,000 yen/ha/yrto 20,000
yen/ha/yr seemed to cover the payment from the rotation
age of 35 to 75 years.

CONCLUSIONSIn this paper, we conducted cost evaluation of the snow

Fig. 4. Difference in NPV.

damage protection through the optimization framework
for the forest stand management of sugi forests. If theTable 3. Annual payment under different rotation age
Rotation Age

Annual Payment(yen/ha/yr)

induced costs by the snow protection management scheme

20

9869

becomes more than the expected cost of snow damage on

25

85900

the forest management, we would stay with the optimal

30

60385

thinning regime without considering snow damage pro-

35

20626

tection into the management scheme. If not, it would be

40

15060

better off for the forest manager to maintain the optimal

45

12102

regime under the snow protection scheme. In the analysis,

50

10416

55

10503

we utilized the height-DBH ratio for the management

60

19075

65

19006

70

17723

75

16758

scheme against snow damage for the cost evaluation ofsnow damage. The previous study indicated that the H/Dratio should be maintained less than or equal to 65%
against snow damage, so that we used this requirement asa constraint into the optimization framework.

rotation length, so that the following annual payment was

rotation age was 65 years for the unrestricted case, whileused to compare the results.
1 T
1− (
)
+
r
1
CSD = APT
1
1− (
)
1+ r

Our experimental analyses showed first that the optimal
50 year rotation age was for the restricted case. The
difference in NPV of these two cases as the induced costs(15)

ranged from 179,867 to 1,910,713 yen/ha over the rotation age of 20 to 75 years. Under the optimal rotation of
65 years, the cost became 914,226 yen/ha. In order towhere
CSD : Induced cost for snow damage protection

reduce an effect of different length of the rotation age,

APT : Annual payment

annual payment was next calculated based on the dif-T

: Rotation age

ference in NPV. The largest was 85,900 yen/ha/yr of ther: Interest rate (1%)

age 25 years, while the smallest was 9,869 yen/ha/yr of
the rotation age 20 years. At the optimal rotation age ofTable 3 shows the resultant annual payment under

50 years for the restricted case, it was 10,416. All in all,different rotation ages. As was observed in Figure 4,

10,000 yen/ha/yr to 20,000 yen/ha/yr seemed to cover theannual payment was the largest for the rotation age of 25

payment from the rotation age of 35 to 75 years.

years. It ranged from 9,869 yen/ha/yr of the rotation age

Although the results indicated here are limited to theEconomic Analysis of Snow Damage on Sugi (Cryptomeria japonica) Forest Stands in Japan Within the … ‧ 149

assumed situations, the results suggested that cost evaluation through the provided approach is useful for the forest
policy toward snow damage protection through the management scheme. The induced costs can be used as thebasis for the coverage of snow damage in the forest management. The cost for forest owners could be offset by thiscoverage.

LITERATURE CITED
Adams, D. M., and E. K. AR. 1974. Optimizing the managementof unevenaged orest stands. Can. J. Forest Res. 4: 274-287.
Amidon, E. L., and G. S. Akin. 1968. Dynamic programming todetermine optimum levels of growing stock. For. Sci. 14: 287291.
Arimizu, T. 1958a. Regulation of the cut by dynamic programming.
J. Oper. Res. Soc. Japan 1: 175-182.
Arimizu, T. 1958b. Working group matrix in dynamic model of forest
management. J. Jpn. For. Soc. 40(5): 185-190.
Arthaud, G. J., and W. D. Klemperer. 1988. Optimizing high andlow thinnings in loblolly pine with dynamic programming. Can.
J. For. Res. 18: 1118-1122.
Brodie, J. D., and R. G. Haight. 1985. Optimization of silviculturalinvestment for several types of stand projection systems. Can. J.
For. Res. 15: 188-191.
Brodie, J. D., and C. Kao. 1979. Optimizing thinning in Douglas-fir
with three descriptor dynamic programming to account for accelerated diameter growth. For. Sci. 25(4): 665-672.
Brodie, J. D., Adams, D. M. and C. Kao. 1978. Analysis of economic
impacts on thinning and rotation for Douglas-fir, using dynamic
programming. For. Sci. 24(4): 513-522.
Chen, C. M., Rose, D. W. and R. A. Leary. 1980. Derivation of
optimal stand density over time-a discrete stage, continuous state
dynamic programming solution. For. Sci. 26(2): 217-227.
Haight RG, Brodie JD, Dahms WG (1985) A dynamic programming algorithm for optimization of lodgepole pine management.
For. Sci. 31(2): 321-330.
Hiller FS, Lieberman GJ (1990) Introduction to operations research.
McGraw-Hill, New York.
Kato A, Nakatani H, Taira H (1992) Relationships of snow damage
to stand and topographic factors in sugi stands. J. Jpn. For. Soc.

74(2): 114-119. (In Japanese).
Kato, A. and K. Tanaka. 1997. Parameter estimation for Boka sugi
forest stands in the system growth simulato. Chubu Forest Research
45: 43-46 (In Japanese).
Kato, A., Zushi, K. and K. Tanaka. 2008. The growth parametersof the system yield table “Silv-no-mori” in Kawaidani-sugi clone(Cryptomeria japonica D.Don). J. Toyama FFPR Res. Center. 21:
9-16 (In Japanese).
Martin, G. L. and E. K. AR. 1981. A dynamic programming analysis
of silvicultural alternatives for red pine plantations in Wisconsin.
Can. J. For. Res. 11: 370-379
Paredes, V. G. L. and J. D. Brodie. 1987. Efficient specificationand solution of the evenaged rotation and thinning problem. For.
Sci. 33(1): 14-29.
Riitters, K., Brodie, J. D. and D. W. Hann. 1982. Dynamic programming for optimization of timber production and grazing in ponderosa pine. For. Sci. 28(3): 517-526.
Schreuder, G. F. 1971. The simultaneous determination of optimal
thinning schedule and rotation for an even-aged forest. For. Sci.
17(3): 333-339.
Tanaka, K. 1995. A system for prediction a forest stand table andheight curve (Silv-Forest). In: Research on development of growthsimulators, Konohira, Y. (ed.). pp. 22-32. (In Japanese)
Torres-Rojo, J. M. and Brodie, J. D. 1990. Demonstration of benefits
from an optimization approach to the economic analysis of naturalpine stands in Central Mexico. Forest Ecol. Manag. 36(2-4): 267278.
Valsta, L. T. and J. D. Brodie. 1985. An economic analysis of hardwood treatment in loblolly pine plantations-a whole rotation dynamic programming approach. In: The 1985 symposium on systemanalysis in forest resources, Dress PE, Field RC (eds.). Georgia
Center for Continuing Education, Athens, pp. 201-214.
Vonesh, E. V. and R. L. Carter. 1987. Efficient inference for randomcoefficient growth curve models with unbalanced data. Biometrics
43: 617-628.
Yoshimoto, A. and R. Marušák. 2007. Evaluation of carbon sequestration and thinning regimes within the optimization frameworkfor forest stand management. Eur. J. For. Res. 126(2): 315-329.
Yoshimoto, A., Paredes, V. G. L. and J. D. Brodie. 1988. Efficient
optimization of an individual tree growth model. In: The 1988
symposium on systems analysis in forest resources, Kent BM,
Davis LS (eds.). USDA Forest Service. General Technical Report
RM-161, pp. 154-162.