Demonstrating the Effect of Height Variation on Stand-Level Optimization with Diameter-Structured Matrix Model

: The weakness of the population matrix models is that they do not take into account the variation inside the class. In this study, we introduce an approach to add height variation of the trees to the diameter-structured matrix models. In this approach, a new sub-model that describes the height growth of the trees is included in the diameter-structured model. We used this height- and diameter-structured matrix model to maximize the net present value (NPV) for the remaining part of the ongoing rotation for Scots pine ( Pinus sylvestris L.) stand and studied how the height variation affects to the results obtained through stand-level optimization. In the optimization, the height variation was taken into account by setting the lower saw-log price for the short trees. The results show that including the height variation into the optimization reduced the ﬁnancial outcome by 16–18% and considerably changed the structure of optimal management (e.g., timings for thinnings, rotation period and intensity of thinnings). We introduced an approach that can be applied to include not only height variation but also variation of other tree properties (such as branchiness or the amount of heartwood and sapwood) into the matrix models.


Introduction
Structured population models (either size-, stage-or age-structured) provide the connection between the population-level dynamics and individual level vital rates [1]. In brief, for stage-or age-structured populations a specific model is constructed to describe reproduction, growth and mortality [2] with the probabilities that population in one age or stage (size) class transfers to another between two sequential periods in time [3,4]. Growth is modeled as a transition from a class to upper classes, survival (opposite to mortality) as the cumulated transitions from a class to other classes and recruitment (reproduction) as a transition into the first class [5]. Particularly in forestry size-structured models (e.g., interchangeably transition matrix model, matrix population model, matrix growth model or matrix model) are thriving with applications covering a wide range of areas [6,7]. For instance, transition matrix models have been applied in discovering biodiversity dynamics [8], in studying the impacts of climate change [9], in biomass carbon stock studies (e.g., [10]) and in assessing optimal forest management [11][12][13][14]. For optimal forest management, there are also other approaches widely applied for predicting tree growth-such as process-based [15][16][17][18] and statistical-empirical growth models [19][20][21]. One argument favoring matrix models over e.g., individual-based statistical-empirical models is the fact that they in most cases require distinctively less amount of information [22].
Despite the extensive application of matrix growth models they have been criticized for their inability to incorporate variation among the individuals consisting of a particular size or age category [23]. Since matrix models operate on size-(or age-, or stage-) structured populations differences among individuals due to size (age or stage) might be an incomplete predictor of actual growth [7]. Further, the outputs of the matrix models are sensitive to the dimension of the matrix, or equivalently to the width of the size classes. This dependence of elasticities on matrix dimensionality has been criticized to lead possible fast pathways which in turn may result in biologically and/or physically impossible growth rates [7]. However, there are methods to circumvent this problem [7]-nevertheless, they are tedious, and require more information which to some extent violate Occam's razor principle to apply parsimonious models [24].
For modeling tree breeding with a matrix model, it is a problem that variation between individuals within the same size (age or stage) class is not taken into account. Namely, a recent study [25] suggests that branch properties (i.e., branch angle, branch thickness and branch multiplicity) are statistically significantly better in improved trees (i.e., genetically improved by tree breeding methods) than in unimproved trees. However, there is considerable variation between individuals. Stated differently, valuable branch traits (flat branch angle, thin branches and fewer branches in a whorl) are more common with improved material than with unimproved, and this can further be converted into monetary value [26]. Currently, more studies focusing on combining branch traits with lumber grading and further pricing [27,28] are called for. Further, more research on branch characteristics associated with unimproved and improved seed material [29,30], and their effect on timber grade distribution is needed. However, before such results exist, one can estimate the effect of fluctuating tree characteristics on financial performance by conducting simulations.
This study introduces an approach that enables us to take into account the height variation between individuals within a particular diameter class in population matrix models. It applies state constrained optimal control problems related to forest management. The state is described with a diameter-structured population matrix model, where is included the new sub-model describing the height growth of the trees. The computational methods of this research are based on the numerical methods of nonlinear, constrained optimization and differential equations, which previously followed by Pyy et al. [14]. The main aim of the study is to demonstrate possible effects of height variation on financial performance when a diameter-structured matrix model is applied in stand-level optimization.

The Optimization Problem
Typically, the optimization helps in solving the forest stand management policy that maximize the net present value (NPV) of the ongoing rotation. In order to formulate the control problem we define the following notations. We divided trees into non-overlapping diameter and height classes with uniform class widths ∆x and ∆z, respectively. The diameter of the trees in the diameter class i = 1, . . . , N is denoted by x i and the height of the trees in the height class j = 1, . . . , K by z j . Time was divided into subintervals [t k−1 , t k ), k = 1, . . . , M, with the length of the time step ∆t. We denoted a tree in the diameter class i and the height class j tree (i, j). Denoting number of trees (i, j) and number of removed trees (i, j) per unit area at time step k by y k ij and h k ij , respectively. We denoted y k = (y k 11 , . . . , y k 1K , y k 21 , . . . , y k 2K , . . . , y k N1 , . . . , y k NK ) and . By y and h are denoted vectors (y 1 , . . . , y M ) and (h 1 , . . . , h M−1 ) , respectively.
The optimization problem is in form where w k ij is the discounted revenues of a tree (i, j) at time step k. The clearcut is done at time step M. We solved the optimization problem with different values of M (time step ∆t remains same so final time t M changes) and chose that which gave the maximal value for the problem (1). Every tree is divided into two quality classes: saw log and pulpwood denoted by subindexes q = 1 and q = 2, respectively. Further, if the height of a tree in diameter class i, is less than µ i − 0.5 m then the unit price of the saw log for the tree is lower than the original saw log price, where µ i denotes the mean height of the diameter class i. Technically, the term µ i − 0.5 m is a sanction which degrades the saw log price for trees which are shorter (by 0.5 m) than the mean height of the particular diameter class i. The rationale is that the smaller the tree the less allowance there is for bucking, which in turn has a negative effect on a saw log price [31]. In this connection it should be highlighted that although the height variation does not have a direct correspondence to wood quality it has bearing on the saw log recovery. The discounted price is where V ijq is the volume of that part of the tree (i, j) that belongs to quality class q, c 1 and c 2 are the price of saw log and pulpwood per volume, respectively, and c 1 is the saw log price for trees shorter than µ i − 0.5 m. Technically, three alternative unit prices for saw logs were applied for trees shorter than µ i − 0.5 m. The applied unit prices were: (a) c 1 = c 1 , (b) c 1 = 0.8c 1 and (c) c 1 = c 2 . Price for saw logs (c 1 ) was 58.44 em −3 and for pulpwood (c 2 ) 16.56 em −3 .
State y k depends on number of removed trees h k via matrix population model [13] where R k is number of new trees at time step k and y 0 is a size distribution of trees at time t 0 . The sets of constraints for y and h are respectively, where h max is the upper limit of removed trees and B = 50 m 3 ha −1 denotes the lower limit of making a thinning at time step k.
The matrix A k is in form The coefficients a k i and b k i are assumed to depend on basal area y k i of the stand. The details about model was described by Pyy et al. [13]. The matrix H is a K × K-matrix that describes how the trees move from a height class to another. We wanted to choose a small height class width ∆z so we assumed that trees from the height class j can remain in the same class or grow to every upper height class j (j ≥ j).

So the matrix H is in form
where 0 ≤ e j j ≤ 1, j , j = 1, . . . , K are proportions of trees that move from the height class j to the height class j . We assumed that H does not depend on time step k or diameter class i. The parameters are fitted on the data described in Section 2.2 below followed by other study [13]. The projected gradient method was used to solve the optimization problem (1) in accordance with other study [14].

Model Estimation and Data
The data used to estimate the parameters for the population model were derived from two long-term experiments (i.e., HARKAS series; for detail see, e.g., [32]). The two experiments were established in the 1978 and 1984 in even-aged, pure commercial Scots pine (Pinus sylvestris L.) stands located on mineral soil in Muhos municipality, Ostrobothnia region, Finland (26 • 06 05 E and 64 • 46 02 N, asl 60-70 m). The data were based on 8 360 tree records from 20 sample plots. The stand management among the sample plots fluctuated considerably: from control plots (no thinnings) to very intensive thinnings (60% of the basal area removed). The stand age in the data set was on average 25 years, indicating that establishment costs and costs of precommercial thinning are sunk. Further, we did not apply any modeling for the time period between stand establishment and 25 years.
We fitted the tree height data (2205 trees, heights between 2.7 m and 19.3 m) to follow Näslunds height curve [33] x Parameters b 0 and b 1 were estimated by least-squares method. We assumed that at time t 0 height of the trees in the diameter class i follows normal distribution [33] with mean and standard deviation where σ 2 is the variance of the variable The minimum diameter of quality class q is denoted by x min,q , q = 1, 2. For saw log x min1 = 15 cm and for pulpwood x min2 = 5.5 cm [28]. Let l ijq denote the height, where the tree (i, j) has the diameter x min,q . Then for the tree (i, j), the saw log part is from ground to the height l ij1 and pulpwood part from height l ij1 to the height l ij2 . The height l ijq = (1 − s)z j , where 0 ≤ s ≤ 1. The unknown s was solved from the equation [34] x min,q f (s j ) In the equation s j = 1 − 1.3 z j and f (s) is taper curve where the form and the values for parameters b 1 , . . . , b 8 were taken from [34]. We solved s from the Equation (13) by using MATLAB function SOLVE.
The volume of a tree (i, j) is where the form and the values for parameters v 1 , . . . , v 5 were taken from [34]. For a tree (i, j) the area under the taper curve (14) on the intervals (0, z j ), (l ij0 , l ij1 ) and (l ij1 , l ij2 ) are A ij0 , A ij1 and A ij2 , respectively, where l ij0 = 0. Then the volume In the population model (3) the parameters a i and b i of the matrix A k were estimated similarly as in explicit approximation of continuous diameter structured population model (see [14]). The parameters e jj , j, j = 1, . . . , K, of matrix H were estimated with the proportion estimator [6]. The recruitment model (4) is based on the average stand conditions in the 7th National Forest Inventory [35]: where P pine = 0.5581 is relative proportion of pines, S re f = 4009.5 is average number of trees per hectare, S = ∑ N i=1 ∑ K j=1 y k ij is actual number of trees per hectare and d is number of time steps before the new trees are big enough to the first size class. Coefficient describes the density effect at time step k as a function of stand basal area BA k . Finally, in the simulations, the time step was set to 5 years with width of diameter class ∆x = 3 cm and width of the height class ∆z = 0.5 m. Minimum diameter and height were set to 2 cm and 2 m, respectively. Then, a minimum of 15-yr growth was assumed to new trees to achieve the first size class, d = 3 cm.

Results
The simulations revealed that height variation combined with pricing had a considerable effect on both financial performance as well as on the optimal management schedules (Table 1, Figure 1). For instance, when saw log unit price was set equal to pulpwood price for trees shorter than µ i − 0.5 m the optimal rotation period was shortened. Especially when discounting with 4%, the optimal rotation period was shortened 20 years ( Figure 1). Further, when discounting with 3%, the structure of thinnings changed drastically with saw logs price set to equal pulpwood price for trees shorter than µ i − 0.5 m, compared to the other two price options ( Table 1). The financial performance (expressed as maximum net present value, MaxNPV) decreased significantly, ranging from 16% to 18% (depending on the interest rate) when height variation with pricing was taken into account ( Table 1). Regardless of the discount rate (3% or 4%) height variation with pricing shortened the optimal rotation period but also reduced the financial outcome. When sawlog unit price was set equal to pulpwood unit price the MaxNPV reduced 409-477 e. However, the impact of height variation and pricing to mean annual increment (associated with stand-level optimum) was minor: approximately 1% difference among the pricing options (Table 1).  As an example of how height variation develops in a diameter class, we will consider the case where the interest rate is 3% and the price of saw log for short trees is the price of saw log.

Discussion
Forest growth and yield models and optimization models are both needed in order to achieve the goals set for forest management [36]. There are, of course, several alternatives for both growth and yield models [37] and optimization models [36,38] to be applied. Currently, traditional empirical-statistical models have made way to e.g., process-based models [36] and to matrix models [6]. Despite the fact that matrix models have been applied to almost all subject areas of forestry [6], they have been criticized for their inability "to incorporate variation among individuals within a size category" [7]. This deficiency associated with matrix models is particularly troublesome regarding tree breeding since it has recently been shown that improved reforestation material tends to have a different distribution with respect to e.g., wood properties (such as branchiness) than unimproved material [25]. Since distributional properties are likely to be different among improved and unimproved reforestation material, this leads to a question of whether that might have an impact on financial performance as well. To study this, the constraint of homogeneity associated with individuals within a particular size category (cf. [2]) needed first to be relaxed. The main objective of this study was to introduce a variation between individuals into a matrix model, and further to demonstrate its effect on financial performance assessed according to stand-level optimization (cf. [39]).
Technically, the height variation between individuals was included in the financial analysis (stand-level optimization) by setting a penalty. In brief, a tree shorter than an arithmetic average tree height within a size category was labeled with a lower unit price of saw logs. The rationale was that the shorter the tree the less allowance there is for bucking, which in turn has a negative effect on a saw log price [31]. Our results indicated that including variation into the matrix model has a notable impact on financial performance. Interestingly, the MaxNPV changed by approximately 16-18% when height variation with pricing was included in stand-level optimization. Although the variation was here set to tree height (and further linked to saw log pricing by specific penalty function), it does not matter which variation is eventually included in the matrix model. One can choose variation in any stem property, for instance branchiness [40] as long as the variation can be linked to financial valuing. Further, the ideal case would be to include variation in some of the key indicators of timber quality-such as the size and type of knots or the amount of heartwood and sapwood [37]-into the matrix model, and to examine this effect in monetary terms with and without the genetic gains [26]. However, such an analysis strongly depends on available data sets, and thus it is conditional to status quo-whether such data had been gathered until now or not.
Finally, it should bear in mind that this study excludes any consideration of sustainability (cf. [41]) or ecosystem services (cf. [10]), solely focusing on the financial performance of timber production. Since this study introduces a new method to include the variation between individuals in a matrix model, and further to assess it in monetary value, the authors consider this contribution to be the first step, and further steps in the near future would cover more than merely timber production.

Conclusions
Including height variation of trees in a diameter-structured matrix model changes the optimal solutions compared to the case without height variation. Differences (with respect to rotation period and intensity as well as the timing of thinnings) are clear, regardless of the applied discount rate. This simply indicates that variations in tree characteristics (e.g., in-branch traits) within the same size class needs to be taken into account more systematically than has been done earlier in population matrix models (e.g., [7]). The approach used in this study enables including any property relevant to timber grading and further pricing and maximizing financial performance.