A Forest Growth Model for the Natural Broadleaved Forests in Northeastern Korea

Abstract

While a large sum of timber stock in private forests, especially broadleaved forests, has been ignored by their owners, a rising global concern about climate change and ecosystems has led to a renewed interest in natural broadleaved forest management strategies. This study establishes the forest growth model for the natural broadleaved forest of Gangwon-do based on the matrix model developed by Buongiorno and Michie. The matrix model by Buongiorno and Michie has been widely applied to study forest population dynamics, especially for uneven-aged forests. To develop an existing matrix model, our approach applies transitional probabilities of forest stands which are calibrated using National Forest Inventory data. Both long and short-term predicted simulation results show that the predicted average tree density and diameter distribution from our model are very close to the stand density and diameter distribution from observed data. Although the model simplifies reality, the results from our study confirm that our models are valid enough to predict the average stand status of the broadleaved forests in Gangwon-do.

1. Introduction

In Korea, about 64% of the total forestland is privately owned. The private forests fail to attract forest owners’ concern due to low productivity and small-scale ownership despite their large sum of accumulated timber stocks. While these forests have been far from the major concern of foresters, recent development in the field of climate change has led to a renewed interest in natural forest management strategies. Thus, it is becoming difficult to neglect the existence of unmanaged private forests, and appropriate tools are needed to support forest management decision making for these forests. Such a necessity becomes even more significant when it comes to deciduous forests in which unmanaged practices are more prevalent than with any other types of forests, but also in considering that an additional 39% of the total forest area is predicted to change into broadleaf forests in a century [1]. This study employs a quantitative model approach to evaluate forest growth dynamics of natural broadleaved forests in Ganwon-do (province). Ganwon-do is located in northeastern Korea, over the latitudes 37°02′ to 38°37′ N and the longitudes 127°05′ to 129°22′ E. This region is a mostly mountainous area with affluent forest resources, and broadleaved forests account for 37% of the total forests in the area. The major climate characteristic of the area is defined as a hot and humid continental climate according to Köppen-Geiger’s classification, but warm, humid climatic conditions and humid subtropical climatic conditions locally appear near the shoreline. A growing body of literature has investigated forest growth dynamics because forest management actions are based on information about current and future forest dynamics (i.e., growth, mortality, reproduction, and related changes in the stand) [2]. From a management perspective, forest growth models are crucial because they provide information for decision making through updating inventories, predicting future yields, and exploring management alternatives [3]. Individual-based models and matrix models are common approaches for studying forest dynamics. Individual-based models can explain the spatial location of individual trees and account for spatial heterogeneity [4]. However, because individual-based models need to deal with a much larger amount of information than matrix models [5], it is impossible to apply them to large-scale studies [6]. A matrix model was initiated to study populations of animals and vegetation communities [7,8] and adopted to study forest populations [9,10]. The matrix model is now one of the most widely used tools for studying forest growth. Matrix models are useful for modeling the size distribution at the population level rather than individual trajectories, especially when the predicted population has many individuals [6]. Another widely applied forest growth model is the GAP model which simulates forest development by summing the dynamic of some regeneration patches [11,12]. However, since the GAP model was developed to study natural developments of tree stands, it is inappropriate to study the growth and yield of managed stands [13]. Bartelink [14] investigated growth dynamics of mixed forest stands in Western Europe using the COMMIX (COMpetition in MIXtures) model. The simulation results show that the COMMIX forest model is fitted to evaluate the development of monospecies stand as well as mixed stands of Douglas-fir and beech. However, the model needs to incorporate other features such as natural regeneration and water related relationships.

Despite increasing interests in the use of the forest growth model, there has been little discussion about the potential use of the dynamic forest growth model for the natural broadleaved forests in Ganwon-do. Moreover, most domestic studies have been far more focused on the even-aged, mostly pure, coniferous forests, because of data availability, a lack of motivation due to relatively low profitability, and problems related to complicated forest structures [15]. Thus, the nature of natural deciduous forests in Gangwon-do remains unclear. Shin [16] examined the environmentally friendly forest management methods of nationally owned natural deciduous forests in Pyungchang, Ganwon-do, and considered factors concerning biodiversity, but he has not dealt with dynamic forest growth but rather used observed static data to evaluate forest growth. The lack of consideration of dynamic features of forests could result in unsound deciduous forest management strategies, which in turn could exacerbate less than ideal forest conditions and accelerate the current use of unmanaged practices.

To comply with current political and academic needs of describing forest dynamics of the private owned natural broadleaved forests in Ganwon-do, we adopt a matrix model developed by Buongiorno and Michie [10] using transition probabilities. This model was designed for uneven-aged forest management with a selection method. The matrix model was applied for studying uneven-aged forests in the Valsugana valley of the Trentino province, and the model provided accurate short-term forecasts and reasonable results in long-term simulations of forest dynamics [17]. Several studies, such as [18,19], developed forest optimization models based on Buongiorno and Michie [10]‘s work, but adding tree size [18] and species diversity variables [19]. Other studies [20,21] include carbon sequestration service to the matrix model for considering climate change issues in forest dynamics. In the general form of matrix model, parameter estimations are diversified by two aspects: (1) state dependent (i.e., non-linear programming); (2) state independent (i.e., linear programming). In general, trees grow faster and die less when there is low competition in the nearby area [22]. Therefore, the growth and survivor rates of trees in the forest are inevitably affected by changing forest density (i.e., basal area and the number of trees per a unit of area). Harvesting can bring the same impacts because removing trees can create more room for trees in the nearby area. Some studies have attempted to incorporate this relation into the model by setting transition parameters as variable rather than fixed [19,23]. However, it is not always necessary for the matrix growth modeling to incorporate state dependency [24] because the observed biological changes resulting from area competition could be too small to be statistically significant [20]. Based on these theoretical backgrounds, we assume the changes in stand density (i.e., basal area, the number of trees per hectare) do not significantly affect the stand growth in that a diameter growth is mainly affected by heights of neighbouring trees rather than basal areas of them in natural broadleaved forests [25,26]. In this study, we predict average stand growth over the province using the matrix model based on transitional probabilities of forest stands. We provide predicted simulation results and their statistical evaluation in the short and long terms. To calculate the transitional probabilities, we used data from the Korean national forest inventory (NFI).

2. Data Preparation for the Model

2.1. National Forest Inventory (NFI)

The data was collected from the Korean NFI. According to data from the NFI, a total of 822 plot-clusters with 2877 plots are located in Gangwon-do. The NFI only offers the overall stand information (e.g., how many trees per a plot are in the first size-class), so, it is hard to identify the individual tree growth (e.g., how much an individual tree grows during θ years).

The sampling design is a 4 km × 4 km grid, and a permanent sampling plot cluster is located at every grid point. Each year, about 20 percent of the entire plots are measured. Therefore, the entire measurements are completed in five years. In general, a cluster consists of four sub-plots (Figure 1), but some clusters are made up of less than four sub-plots (hereafter “plot”). The plot is composed of four circular zones with different scales for each subject matter; infant trees, small and medium trees, large trees, soil, and other vegetation. The trees in small and middle DBH (Diameter at Breast Height) classes (6–30 cm) are measured within an 11.3 m radius (0.04 ha) and those in large DBH classes (30 cm~) within a 16 m radius (0.08 ha). Meanwhile, dead trees and logged trees are only measured in the middle circular zone of an 11.3 m radius regardless of their tree sizes and are reported only from the center of the cluster (hereafter “S1 plot”). Consequently, there is less information available for the dead trees.

2.2. Data Collection

In line with the objective, the data were collected from (1) unmanaged and (2) broadleaved forests within the period of 2006–2009 and 2011–2014. The unmanaged forest is defined as the forest that has naturally regenerated without planting and has not been logged for last 10 years. The broadleaved forest is defined as the forest where broadleaf species are predominant, representing more than 75% of canopy cover or tree density. Initially, 1106 plots classified as ‘natural’ and ‘broadleaved’ forests were collected. After some filtering process (see the below list), it reduces the total to 134 plot clusters with 441 plots which are eventually used to estimate transition probabilities, $a_{i}and{b}_i$. The number of living trees and their DBHs were obtained from the remaining plots.

Below is the list of the filtering criteria:

• the plot clusters where logging practices for the last 10 years are reported in the S1 plot
• the plots that have missing data of the fifth NFI
• the plots that tree distribution is incorrectly reported
• the plot clusters with only one plot

The list of dead trees is separately recorded from that of living trees since they are only measured in the S1 plots. For the mortality rate estimation, the similar but slightly different filtering criteria were applied due to its small size:

• the plots where logging practices for the last 10 years are reported in the S1 plot
• the size classes that have more dead trees than living trees

On the other hand, the missing data of the fifth NFI is recovered (only for S1 plots) to ensure enough sample size for the estimation. As a result, a total of 132 S1 plots remain to estimate the mortality rate. Although almost half of the plots are filtered out, the remaining plots are evenly distributed over the province compared to the original plot distribution (see Appendix A). Therefore, the random distribution assumption of plots still holds. The plot summary of 441 plots is shown in

Table 1. Summary of the 134 clusters with 441 subplots of the sixth national forest inventory (NFI).

	Variable	Mean						Range		SD
Plot characteristics	Stand ages (years)	48.89						21~103		11.1
	Stand density (tree/ha)	928.82						138~2444		404
	Basal area (m2/ha)	23.58						2.3~51.5		10.3
	Crown density Distribution (%)	0%–25%			26%–50%			51%–75%		76%–100%
		1			6.3			25		68
Topography	Elevation (m)	728.57						45.0–1383.0		280.4
	Slope (°)	34.25						5.0–60.0		9.2
	Aspect Distribution (%)	E	N	NE		NW	S	SE	SW	W
		17	14	15		10	12	10	9	12

.

3. Forest Growth Model

3.1. A Matrix Model Structure

In the matrix model, the stand state is described by the number of trees in the finite number of diameter classes (i = 1, n). The column vector denotes the number of living trees, per a unit of area, in diameter class i, at specific time t as follows:

$${\underset{\_}{\mathit{y}}}_{\mathit{t}}=\left[y_{it}\right].$$

(1)

The number of harvested trees, per a unit of area, in diameter class i at specific time t during the period $\Theta ,$ is denoted by the column vector:

$${\underset{\_}{\mathit{h}}}_{\mathit{t}}=\left[h_{it}\right].$$

(2)

Subsequently, the status of forest stand at time $t+\Theta$ can be presented as the following n equation systems (3):

$$\begin{array}{c}{\mathit{y}}_{\mathbf{1},\mathit{t}+\mathit{\theta }}=\mathit{I}+{\mathit{a}}_{\mathbf{1}}\left({\mathit{y}}_{\mathbf{1}\mathit{t}}-{\mathit{h}}_{\mathbf{1}\mathit{t}}\right),\\ {\mathit{y}}_{\mathbf{2},\mathit{t}+\mathit{\theta }}=b_2\left({\mathit{y}}_{\mathbf{1}\mathit{t}}-{\mathit{h}}_{\mathbf{1}\mathit{t}}\right)+a_2\left({\mathit{y}}_{\mathbf{2}\mathit{t}}-{\mathit{h}}_{\mathbf{2}\mathit{t}}\right),\\ \dots \dots \dots \mathrm{..}\\ {\mathit{y}}_{\mathit{n},\mathit{t}+\mathit{\theta }}=b_n\left({\mathit{y}}_{\mathit{n}-\mathbf{1}\mathit{t}}-{\mathit{h}}_{\mathit{n}-\mathbf{1}\mathit{t}}\right)+a_n\left({\mathit{y}}_{\mathit{n}\mathit{t}}-{\mathit{h}}_{\mathit{n}\mathit{t}}\right).\end{array}$$

(3)

where $a_i$ is a probability that a living tree at time t stays alive at time $t+\Theta$ in the same diameter class. The $b_i$ is a probability that a living tree in a diameter class i–1 at time t advances to a larger diameter class i at time $t+\Theta$. Lastly, the vector $I$ represents the ingrowth function, which is the number of expected trees entering the smallest diameter class during the growth period $\Theta$.

The ingrowth function is defined as follows:

$${\mathit{I}}_{\mathit{t}}={\beta }_0+{\beta }_1{{\displaystyle \sum }}_{i=1}^n{\mathit{B}}_i\left({\mathit{y}}_{\mathit{i}\mathit{t}}-{\mathit{h}}_{\mathit{i}\mathit{t}}\right)+{\beta }_2{{\displaystyle \sum }}_{i=1}^n\left({\mathit{y}}_{\mathit{i}\mathit{t}}-{\mathit{h}}_{\mathit{i}\mathit{t}}\right).$$

(4)

where the sign of each parameter is expected to be: ${\mathit{I}}_{\mathit{t}}\ge 0,{\beta }_1<0,{\beta }_2>0,{\beta }_0\ge 0$. The ${\mathit{B}}_{\mathit{i}}$ is the basal area of the tree of an average diameter of each diameter class i. The parameter ${\beta }_1$ represents the effect of the total basal area of the stand on ingrowth. This parameter is expected to be negative because dense forest causes higher competition between trees. The parameter ${\beta }_2$ represents the effect of the number of trees of that stand on ingrowth with other things being equal. One might think more trees cause higher stand density and that hence, the parameter should be a negative value. However, having more trees in the stand causes the higher chance of seed dispersal as well. Thus, the larger number of trees will cause a higher rate of ingrowth for a given basal area. The parameter ${\beta }_0$ indicates the possible ingrowth that occurs independently on the stand state, possibly due to seed dispersal from surrounding stands; thereby it is expected to be non-negative.

By replacing ingrowth vector ${\mathit{I}}_{\mathit{t}}$ in Equation system (3) with ingrowth function Equation (4), the new first size class growth function is derived as follows:

$${\mathit{y}}_{\mathit{i},\mathit{t}+\mathit{\theta }}={\beta }_0+d_1\left({\mathit{y}}_{\mathit{i}\mathbf{1}\mathit{t}}-{\mathit{h}}_{\mathit{i}\mathbf{1}\mathit{t}}\right)+\cdots +d_n\left({\mathit{y}}_{\mathit{n}\mathit{t}}-{\mathit{h}}_{\mathit{n}\mathit{t}}\right).$$

(5)

where $d_1=a_1+{\beta }_1{\mathit{B}}_1+{\beta }_2$; $d_i={\beta }_1{\mathit{B}}_{\mathit{i}}+{\beta }_2$, For i > 1.

From Equations (3) and (5), the final form of the matrix model is defined as follows:

$${\underset{\_}{\mathit{y}}}_{\mathit{t}+\mathit{\theta }}=\mathit{G}\left({\underset{\_}{\mathit{y}}}_{\mathit{t}}-{\underset{\_}{\mathit{h}}}_{\mathit{t}}\right)+\mathbf{c}.$$

(6)

where $\mathit{G}=\left[\begin{array}{cccc}{d}_1& d_2& \cdots & d_n\\ b_2& a_2& 0& 0\\ 0& \ddots & \ddots & 0\\ 0& 0& b_n& a_n\end{array}\right],c=\left[\begin{array}{c}{\beta }_0\\ 0\\ ⋮\\ 0\end{array}\right],{\mathit{y}}_{\mathit{t}}=\left[\begin{array}{c}{y}_t\\ y_t\\ ⋮\\ y_t\end{array}\right],{\mathit{h}}_t=\left[\begin{array}{c}{h}_1\\ h_t\\ ⋮\\ h_t\end{array}\right]$, where the matrix G contains the transition probabilities of both up-growth and ingrowth. The constant vector $c$ includes the ingrowth constant ${\beta }_0$.

3.2. Model Calibration and Short-Term Validation

As described earlier, the model parameters were estimated using the fifth and sixth NFI data. A total of 411 plots and 132 S1 plots were used to estimate transition probabilities ($a_{ij},b_{ij},m_{ij})$.

The trees were grouped by diameter size measured at breast height. There are eight size classes with a 5 cm interval, ranging from 6 cm to 40+ cm. Trees between 6 cm–10 cm were grouped in the first size class while all trees above 40 cm were grouped in the eighth size class. Such a classification assumes that no tree grows more than 5 cm during five years based on empirical data that the average annual growth of broadleaf species is 3.03 mm, with the maximum annual increase of 5.06 mm in a second age class [27]. These categories also correspond to the Sustainable forest management guideline [28].

On the other hand, only one group represents all species. According to NFI data of the study area, the six oak species (Q. mongolia, Q. variables, Q. serrata, Q. aliena, Q. acutissima, Q. dentate) account for 50% of the total number of trees and 58% of the total basal area. Another half of the forestlands are occupied by other broadleaf species of more than 100 species. Other broadleaf species account for 47% of the total number of trees and 38% of the basal area. For coniferous species, they have significant differences from broadleaf species in regards to timber prices and carbon emission factors. However, they are not separately grouped due to their minor occupation in the stand. The coniferous species, mainly Pinus densiflora, represent only about 3% of the total number of trees and 6% of the total basal area. Moreover, their occupation is expected to be even smaller in the future because of their shade intolerance traits. That means their representation at steady-state will be even more negligible than the current state [29].

Table 2. Changes in diameter distribution between two inventories.

Diameter Class (cm)		With 441 Plots
		Fifth NFI	Sixth NFI (n/ha)
1	6–10	382.83	346.10
2	11–15	232.21	217.29
3	16–20	143.64	145.80
4	21–25	89.32	97.57
5	26–30	49.36	56.14
6	31–35	27.30	31.28
7	36–40	13.91	15.87
8	40+	16.18	18.77
Total (n/ha)		954.75	928.82
Basal area (m2/ha)		22.13	23.58

shows the forest structure at the time of the fifth and sixth NFI following the classification.

Since the individual tree growth values are not available in the NFI data, the number of trees that are dead, staying in the same size class, and moving to the upper class, were identified by backward calculations. That is, the mortality rate is firstly computed, and then used to derive the transitional probabilities $a_{i}and{b}_i$. The estimation process is further elaborated in Appendix B. Each probability is defined as follows, using a simple proportion.

For the Mortality rate:

$$m_i=\frac{\mathrm{The}\mathrm{number}\mathrm{of}\mathrm{dead}\mathrm{trees}\mathrm{of}\mathrm{diameter}\mathrm{class}i\mathrm{during}5\mathrm{year}}{\mathrm{The}\mathrm{number}\mathrm{of}\mathrm{living}\mathrm{trees}\mathrm{of}\mathrm{dimater}\mathrm{class}i\mathrm{at}\mathrm{the}\mathrm{fifth}\mathrm{NFI}}$$

(7)

For the probability that a living tree in size class i stays in the same size class during the interval is:

$$a_i=\frac{\mathrm{The}\mathrm{number}\mathrm{of}\mathrm{staying}\mathrm{trees}\mathrm{in}\mathrm{size}\mathrm{class}i\mathrm{during}\mathrm{interval}}{\mathrm{The}\mathrm{number}\mathrm{of}\mathrm{living}\mathrm{trees}\mathrm{of}\mathrm{diameter}\mathrm{class}i\mathrm{at}\mathrm{the}\mathrm{fifth}\mathrm{NFI}}$$

(8)

For the probability that a living tree in i–1 class moves up to the upper class i during the interval is:

$$b_i=\frac{\mathrm{The}\mathrm{number}\mathrm{of}\mathrm{advanced}\mathrm{trees}\mathrm{to}\mathrm{upper}\mathrm{size}\mathrm{class}i}{\mathrm{The}\mathrm{number}\mathrm{of}\mathrm{living}\mathrm{trees}\mathrm{of}\mathrm{size}\mathrm{class}\left(i–1\right)\mathrm{at}\mathrm{the}\mathrm{fifth}\mathrm{NFI}}$$

(9)

The calibration results, including standard error, are provided in Appendix B, but the estimated probabilities are found in the matrix G below.

We estimate the coefficients of ingrowth function using linear regression approach. The dependent variable is the number of ingrowths per cluster, which was obtained from the backward calculation described above. The adjusted coefficient of determination (R²) is 0.242 which means the model only explains a small part of ingrowth variation. Nonetheless, the coefficients are statistically significant at the 0.05 level. More detailed statistical results are found in Appendix C.

The ingrowth function with the coefficients is stated as:

$${\mathit{I}}_t=262.02-10.42{{\displaystyle \sum }}_{i=1}^n{\mathit{B}}_{\mathit{i}}\left({\mathit{y}}_{\mathit{i}\mathit{t}}-{\mathit{h}}_{\mathit{i}\mathit{t}}\right)+0.23{{\displaystyle \sum }}_{i=1}^n\left({\mathit{y}}_{\mathit{i}\mathit{t}}-h_{it}\right)$$

(10)

Finally, the transition matrix G and the constant vector c with coefficients are presented as follows: $\hat{\mathit{G}}=[\begin{array}{cccccccc}0.49& 0.09& -0.03& -0.20& -0.41& -0.66& -0.95& -1.58\\ 0.37& 0.37& & & & & & \\ & 0.33& 0.50& & & & & \\ & & 0.32& 0.54& & & & \\ & & & 0.31& 0.58& & & \\ & & & & 0.31& 0.68& & \\ & & & & & 0.22& 0.69& \\ & & & & & & 0.28& 0.92\end{array}],\hat{\mathit{c}}=\left[\begin{array}{c}262.02\\ \\ \\ \\ \\ \\ \\ \end{array}\right]$.

The short-term predictions, based on estimated probabilities, are compared with the sixth NFI data in

Table 3. Short-term validation.

Diameter Class (cm)		134 Clusters		224 Clusters (No Filtering)
		Observed	Predicted	5th NFI	Observed	Predicted
	Ingrowth	257.54 *	251			260.50
1	6–10	346.10	369.68	435.75	401.53	395.59
2	11–15	217.29	227.56	270.42	257.73	261.28
3	16–20	145.80	148.45	168.39	173.97	173.43
4	21–25	97.57	94.20	98.99	113.45	107.34
5	26–30	56.14	56.32	52.62	61.47	61.21
6	31–35	31.28	33.87	26.25	31.01	34.16
7	36–40	15.87	15.60	13.51	15.49	15.10
8	40+	18.77	18.78	16.37	18.49	18.84
Total (n/ha)		928.82	964.46	1082.30	1073.13	1066.95
Basal area (m2/ha)		23.58	23.96	24.04	25.98	26.00

* Ingrowth is the average number of trees that enter the smallest size group in 441 subplots.

. The forecast slightly over-estimates the number of trees of the small size classes (1 ~ 2), which represents about 5% of the population. Meanwhile, the predictions of other classes (>2) are close to the actual number of trees in each of the size classes. Consequently, the stand density (n/ha) is slightly overestimated due to the differences in small classes while the total basal area per ha is almost accurately predicted. Additional short-term validation was implemented with all non-harvested plots, 549 plots of 224 clusters. The stand density and diameter distribution from extended prediction results are fitted the stand density and diameter distribution from sixth NFI data (The average forest stand of 549 plots shows a higher number of small trees but a lower number of large trees compared to that of 441 plots. As far as the study objective is concerned, this initial stand is not a problem because the steady state is dependent on transition probabilities not the initial state.). Such results confirm that the transition probabilities are valid enough to predict the average stand statues of the broadleaved forests in Gangwon-do.

3.3. Stand Dynamics, Steady-State, and Long-Term Validation

The current forest structure presents a reverse-J shape, which is a typical feature of uneven-aged forests at equilibrium. However, the diameter distribution seems not to reach the steady-state yet. The current state (i.e., the fifth NFI) shows that the large trees (>40 cm) only occupy 13% of the total basal area and 2% of the total number of trees. These proportions, especially the basal area, are too low compared to those of the small classes and ingrowth trees. As large trees dominate the forest stand, their canopies block the sunlight from reaching the understory. Subsequently, it slows down the overall stand growth as the large mature trees substantially restrict the growth of vigorous young trees that are represented here as ingrowth and small-class trees. This process eventually leads to zero forest growth, net of mortality, and the changes of forest distribution would stop at some point which we called a steady-state (or equilibrium). The transition matrix could simulate this process toward a steady-state.

From the growth model, the forest growth during Ɵ is expressed as follows:

$${\underset{\_}{y}}_{t+\theta }=\mathit{G}{\underset{\_}{y}}_t+\underset{\_}{\mathbf{c}}$$

(11)

Based on the assumption that the transition probabilities remain stable over the whole projection period, the forest stand after k Ɵ years from now are defined by the following equation:

$${\underset{\_}{y}}_{t+k\theta }={\mathit{G}}^{\mathit{k}}{\underset{\_}{y}}_0+{\displaystyle \sum }_{i=0}^{k-1}{\mathit{G}}^{\mathit{i}}\underset{\_}{\mathbf{c}}$$

(12)

The stand development over 200 years is presented in (Figure 2). It shows that the number of large trees slightly increases while the number of small trees largely decreases. Figure 3 displays the changes of tree density and the basal area along with stand development. The stand density decreases as the basal area expands, as the basal area directly controls the ingrowth. Moreover, the fluctuation range converges to the equilibrium.

The steady-state requires no change in forest structure during the interval. This state is found at the point where the Equation (12) satisfies the following conditions as indicated in Equation (13):

$${\underset{\_}{y}}_{t+\theta }={\underset{\_}{y}}_t={\underset{\_}{y}}_t\ast$$

(13)

where ${\underset{\_}{y}}_t\ast$ indicates the number of the trees by size class per ha at steady-state.

Table 4. Forest structure at the undisturbed steady-state.

Diameter Class (cm)		Tree Density		Basal Area
		n/ha	%	m2/ha	%
	Ingrowth	111.09
1	6–10	161.92	31	0.81	3
2	11–15	95.10	18	1.26	5
3	16–20	62.76	12	1.60	6
4	21–25	43.66	8	1.81	7
5	26–30	32.23	6	1.98	8
6	31–35	31.22	6	2.67	10
7	36–40	22.16	4	2.51	10
8	40+	77.54	15	13.45	52
Total		526.58	100	26.11	100

below describes the forest structure of the undisturbed steady-state. The most noticeable feature is that the large trees (>40 cm) considerably increase about four times the observed numbers and comprise half of the total basal area. Meanwhile, the ingrowth and small-class trees decrease almost half of the present numbers. The trees between 20 cm and 40 cm diameter change to a small extent. As a result, the size distribution over the diameter classes becomes more balanced than the present state in terms of the number of trees. In addition, it corresponds to the general feature of old-growth forests where the large trees are dominant and oppress the understory growth.

3.4. Long-Term Validation

The predicted steady-state of the broadleaved forest was checked with the old-growth deciduous forest in Kwangnueng, South Korea. The Kwangneung Experimental Forest (KEF) has been designated as a royal tomb forest since 1468, thereby human activities had been completely prohibited until 1913. Under Japanese colonial government, some parts of KEF turned into the experimental forests that exist to date, and the remaining parts have been well conserved [30]. In KEF, the old-growth broadleaved forest is located on Mt. Soribong in the Kwangneung National Forest Reserve that has been protected from human disturbances. Therefore, this forest is considered to be a good example of an undisturbed steady-state of temperate deciduous forests. This forest zone represents cool temperate forests with an annual precipitation of 136.5 cm and temperature of 11.3 °C. The forest species composition and structures were studied in the 100 m × 100 m permanent plot by Lim et al. [29].

Table 5. The Kwangneung national forest permanent plot.

Species	Basal Area (m2/ha)	DBH (cm)		Density (Trees/ha)
		Mean ± SD (cm)	Max.
Quercus serrata	14.16	48 ± 17	110	70
Carpinus laxiflora	6.54	18 ± 15	52	153
Carpinus cordata	2.20	10 ± 7	54	176
Acer mono	1.10	30 ± 17	59	12
Fraxinus rhynchophylla	0.75	33 ± 9	50	8
Sorbus alnifolia	0.56	10 ± 13	50	27
Euonymus oxyphyllus	0.32	3 ± 1	16	530
Cornus controversa	0.25	6 ± 6	38	39
Cornus kousa	0.19	6 ± 4	17	45
Styrax japonica	0.18	4 ± 2	12	100
Acer pseudo-sieboldianum	0.16	4 ± 3	13	89
Styrax obassia	0.11	5 ± 3	15	36
Celtis jessoensis	0.10	11 ± 14	35	4
Prunus mandshurica var.glabra	0.08	9 ± 1	31	6
Others (19 species)	1.32			179
Total	28.00	26.792		1474

presents the old-growth forest stand with tree species and their basal area, mean DBH, and density. In addition, the standing trees are counted from the 2 cm in diameter. While bearing these differences in mind, it is possible to see the general structure of the old-growth forest and compare them with the predicted steady-state, after some counting process. First, four species with mean DBH less than 6 cm are removed from the stand data, which diminishes the density from 1474 trees/ha to 719 trees/ha. Then, half of the number of three species, C. cordata, C. controversa, C. kousa, which have low mean DBH and a wide range of deviation are left out. Such a process is expected to compensate for the loss of some species with the gain of other species. This process gives the stand density of 589 trees/ha and the basal area of 25.91 m²/ha, which is very similar to the predicted steady-state values (i.e., 526 trees/ha and 26.11 m²/ha).

In addition, the old-growth forest shows a similar diameter distribution to that of the predicted forest. The predicted forest suggests that large trees (>40 cm) dominate the forest stands, with a basal area of 13.45 m²/ha and a density of 77.54 trees/ha. Similarly, Q. serrata, the only species that exceeds the mean of 40 cm, accounts for a basal area of 14.16 m²/ha and 70 trees/ha in the old-growth forest. Furthermore, the small number of trees between 20–40 cm are observed in both forests. In the predicted steady-state, the number of trees of each of the size classes between 4 and 7 is less than 10% of total tree numbers, which is noticeably lower than other classes. In the old-growth forest, only two species are found in the range from 20 cm to 40 cm for the average diameter. Given the high deviation, the actual number will be greater than the 20 trees of the two species. However, such results would not change the fact that the old-growth forest tends to have relatively small number of the middle-class trees due to the adverse effects of competition with large trees.

4. Discussion

Our growth model is designed based on the matrix model developed by Buongiorno and Michie. The matrix model has several advantages because it relies on simple linear programming and data requirements: the number of trees in each diameter groups for at least two-time points, thus allowing for aggregate tree data not only an individual format. Therefore, it is relatively easy to calibrate and relatively less dependent on the precision of data. However, these advantages can be turned into a limitation of the model because the model averages out the impact of various site conditions on stand growth. In the real world, the microclimate significantly affects tree growth, especially in mountainous topography [31,32]. Therefore, foresters need to avoid generalizing these approaches, and should instead consider stand characteristics and ecological conditions when they apply the approach to practice. Furthermore, the fixed parameter assumption, which implies that tree growth is not affected by stand structure and basal area, may underestimate the effect of silvicultural management practices, such as thinning and fertilizing, which can promote stand growth [33,34]. One possible alternative is to use a non-linear growth model that allows the dynamic transition parameters of forest evolution [19,23]. However, such an option should be considered when strong correlations between stand condition (e.g., basal area and density) and stand growth are detected [20], because of its complicated calibration process (e.g., appearance of local optima).

Finally, it is recommended that multi-species are incorporated into the model, while our model aggregates all species into one species group. The current model is designed primarily for a subsequent study of optimization strategies of broadleaved forests. Therefore, the major considerations for the classification were (1) timber price and (2) carbon sequestration. Additionally, there were little benefits of classifying species in multi-classes given the facts that (1) the Korean domestic timber market has one single price value for all broadleaved species because most of the broadleaved trees are used as pulp; (2) only six carbon coefficients for broadleaved species have been verified by Korea Forest Research Institute [35,36]. If we were to use the same coefficients to estimate timber volume, carbon stock, and the monetary value of the timber, there would be no benefits of classification considering optimization. Nevertheless, such simplification could cause a significant bias in the analysis in estimating timber and carbon stocks. Also, it could overlook the importance of the diversity of tree species in broadleaved forests management. When more data is available for different species, the current classification should be more highly specified so that we can better understand dynamics of forest growth, especially changes in species composition and their implications on the forest ecosystems.

5. Conclusions

The study develops the forest growth model for the natural broadleaved forest of Gangwon-do. The data collected from the NFI provides sufficient information to estimate transitional probabilities, but its measurement errors (i.e., the number of trees in each size class) were rather high, so negative transitional probabilities at stand level occurred. However, the problem was partly solved by excluding the highly imprecise stands from the estimation. By comparing the predicted average stand growth from our growth model and observed data from NFI, we show that our growth model fits the data well, in both the short term and long term. Applying individual tree growth data (c.f. stand growth data) might be possible alternatives for improving model accuracy, but such individual tree data is currently not available in the Korean NFI. Also, it is recommended to develop multi-species growth model to better address the forest ecosystem and diversity issues.