Survival analyses of individual tree populations in natural forest stands to evaluate the maturity of forest stands: A case study ofpreserved forests in Northern Japan

Abstract

Most survival analyses in forest science have been applied to planted forests where tree age usually had no dispersal on one stand, while there have been few reports about the age-based survival analyses in natural forests. The purpose of this study was to perform survival analyses of individual tree populations in natural forest stands to evaluate the matured states of forest stands. We used a combination of tree-ring and census data from the three preserved permanent plots in pan-mixed and sub-boreal natural forests, Hokkaido, northern Japan. All the living trees (diameter at breast height ≥ 5 cm in 1989) were targeted to identify tree ages using a RESISTOGRAPH. Periodical tree age data with a 10-year age class were used during the observation periods of 1989–1999, 1999–2009, and 2009–2019, and all the changes (i.e., death and new ingrowth) during the periods were recorded. In the analyses, first, we applied survival analyses to find out multi-temporal age distributions and non-parametric estimates. Secondly, we applied parametric Weibull distributions to calculate age-related metrics such as mean lifetime and mean age of stands. Finally, we evaluated these calculated metrics from the viewpoints of matured states of forest stands. The results of non-parametric survival analyses showed the age distribution of multi-modal and exponential shapes. We also found the differences in survival probabilities among periods were not significant except for one plot. We also compared the relationship among estimated mean lifetime, biological lifetime, and mean age of stands derived from parametric survival analyses with the Weibull distribution and evaluated the matured states of stands considering all these aspects. It was implied that the study stands might not get enough matured yet, but some plots showed further progress toward the matured state than others.

INTRODUCTION

A survival analysis, time-to-event analysis, refers to a set of methods for analyzing the length of time until the occurrence of a well-defined endpoint of the event. It is common to use survival analysis in the fields of biology, medicine, engineering, marketing, social sciences, and behavioral sciences (Jung et al., 2012; Gross et al., 2014; Alvi et al., 2015; Ancarani et al., 2015; Buckley et al., 2016). In survival analysis, the covariates of "age" is one of the major parameters to analyze the mortality of target populations accurately (Kalblfeisch and Prentice, 1980; Kleinbaum and Klein, 2011). So, in the field of forest science, most of the survival analyses have been applied to planted forests, where tree age usually had no dispersal on one stand, i.e., no need to identify tree age. In the survival analysis of planted forests, the event is the clear-cut of the stands. Based on this endpoint, Suzuki (1959) has developed his original method to predict the harvesting behavior of planted forests called the Gentan probability theory in Japan. Some other survival analyses had also been applied to planted forests since then worldwide (Waters, 1969; Morse and Kulman, 1984; Amateis et al., 1997; Wyckoff and Clark, 2000).

When it comes to natural forests, the event is the mortality of individual trees. Woodall et al. (2005) carried out survival analyses in natural forests based on DBH instead of tree age. It is generally troublesome to detect tree ages of individual trees in natural forests, particularly in large-scale inventories (Flewelling and Monserud, 2002; Nothdurft, 2013), so most of the researchers used stand metrics such as DBH, basal area, and height to analyze the mortality of trees not only in natural forests but also planted forests (Monserud, 1976; Tatsuhara and Minowa, 1988; Vanclay, 1995; Eid and Tuhus, 2001). In spite of the difficulty in tree age detection, one of the clear advantages in survival analyses with covariates of age, i.e., age-based survival analyses, is to analyze and evaluate the maturity of cohorts considering age distributions, mean age, mean lifetime, etc. In such circumstances, Hiroshima (2014) applied age-based survival analyses to a secondary natural forest in Japan. His study was the first case in natural forests to introduce a survival analysis based on the age of individual trees. He carried out the non-parametric and parametric analyses to illustrate the mortalities of the natural forest stand and tried to evaluate the maturity of stand based on the age-related metrics such as mean age and mean lifetime derived from the non-parametric and parametric estimates. This research was, however, based on only one plot and one-period measurement, so that it was not easy to make a clear general conclusion based on this case study since obtained survival probabilities could be site-specific. Besides, the targeted stand seemed immature, and the estimated survival probabilities seemed not stable over time. Given the current lack of age-based survival analysis, it is required to carry out more case studies of survival analyses, particularly in natural forests. Moreover, it is desirable to apply multi-temporal survival analyses to evaluate the matured states of stands considering the time-series changes in the above age-related metrics.

In this study, we took advantage of the long-term research site in a preserved area located in the University of Tokyo Hokkaido Forest (UTHF). The permanent plots of these sub-boreal forests were ideal for performing the tree age and mortality investigations in the long term because they had not been undergone human-induced disturbances for a long time.

Overall, the purpose of this study was to perform survival analyses of individual tree populations in natural forest stands within the pan-mixed and sub-boreal forests to evaluate the matured states of forest stands. First, we applied survival analysis techniques to study plots to find out multi-temporal age distributions and non-parametric estimates. Secondly, we applied parametric Weibull distributions to calculate age-related metrics such as mean lifetime and mean age of stands. Finally, we evaluated these calculated metrics from the viewpoints of matured states of forest stands. This study was the second case of applying age-based survival analyses for individual tree populations in natural forests in Japan.

MATERIALS AND METHODS

Study Site

The investigation was conducted in the permanent plots of the UTHF (43° 10'–21" N, 142° 23'–41" E), pan-mixed and sub-boreal forests in Hokkaido, Japan (The University of Tokyo Hokkaido Forest, 2019). This forest locates at the transition area from deciduous forests in the cool-temperate zone to coniferous forests in the sub-boreal zone (Kato, 1952) with a mean annual temperature and precipitation of 6.4 °C and 1,297 mm, respectively (The University of Tokyo Hokkaido Forest, 2019). The Stand-based Silvicultural Management System (SSMS), which is a unique natural forest management system conducting the single tree selection of over-matured and defective trees based on natural regenerations, has been practiced since 1958 (Takahashi 2001). Harvesting and management decisions are taken based on the long-term observation data in permanent plots. Among these permanent plots, a total of 25 plots locates in the preserved area, so we call them preserved permanent plots for convenience. These plot sizes range from 0.04 to 2.25 ha with an elevation range between 380–1290 m. Within these plots, DBH measurements of all trees with DBH ≥ 5 cm are performed by UTHF staff regularly, in most cases, 5-years interval. Other than these periodical measurements, basically, there is no human intervention in this preserved area for several decades.

This study chose three samples of preserved permanent plots #5203, 5225, and 5240 within elevation ranges around 600–700 m (Fig. 1 and Table 1). This middle elevation represents the range of similar species composition as well as stand structure with proximity to each other plot in the preserved area. Typical vegetation types in these plots are coniferous and broad-leaved mixed forest dominated by Abies sachalinensis, Picea jezoensis and Acer species (mainly Acer urkrunduense). The main soil types found in these plots were brown forest and podzolic soils (The University of Tokyo Hokkaido Forest, 2019). In addition, these three plots are classified as the same stand type of "coniferous selective cutting with poor regeneration" in UTHF, where we cannot expect sufficient new ingrowth trees continuously. All three plots also have similar slope aspects and slope angles (Table 1) though abiotic factors such as microtopography affected the mortality of trees in the three plots differently. Thus, the selected three plots facilitate the consistent analyses of the study objectives mentioned above.

Fig. 1.

Maps showing the locations of study area and sample plots. (a) UTHF of the Hokkaido island of Japan, (b) preserved area of the UTHF, (c) Contour map of the sampled three preserved permanent plots of UTHF.

Table 1 Stand characteristics of selected plots

Plot #	Plot size (ha)	Elevation (m)	No. of target trees in period 3	Slope aspect	Mean slope angle (degrees)
5203	0.40	580	327	Southwest	15
5225	0.25	690	237	Southwest	20
5240	0.25	600	214	Southwest	25

Collection of Tree Age Data

We used the tree census data (species, DBH, state of alive or dead, cause of death, etc.) of the plots measured between 1989 and 2019. In line with the census data, we set three observation periods of 1989–1999 (period 1), 1999–2009 (period 2), and 2009–2019 (period 3). All the target trees were alive and had DBH ≥ 5 cm in 1989, though parts of these trees were already dead by 2019. We additionally collected the data of tree age after ingrowth (i.e., the number of annual rings from bark to pith in radius at breast height - 2.5 cm, if we set DBH = 5 cm as an ingrowth border) for every target tree. By considering the "age after ingrowth", we can ignore the impact of trees with DBH < 5 cm in the later survival analysis (Hiroshima, 2014). A semi non-destructive device called the RESISTOGRAPH was used to detect the annual rings of both the living and dead target trees (Hiroshima, 2014). The measurement data were then extracted via the DECOM software. Figure 2 depicts the annual rings detected in one of the alive Abies trees after deducting 2.5 cm from the pith using the DECOM^TM software. It can be used to mark tree ring limits automatically or manually based on the boundaries between early and latewood areas (Rinn et al., 1996). RESISTOGRAPH measurements were conducted for a total of 327, 237, and 214 trees in the plots 5203, 5225, and 5240 respectively in 2019 (the end of period 3), and the ages of target trees in period 3 were first identified, and then the ages in period 2 and 1 were calculated by deducting 10 and 20 from those in period 3.

Fig. 2.

DECOM^TM profile showing the drilling resistance of the annual wood growth rings. The differences in density between early and latewood were used to draw the vertical lines.

It was not possible to apply the RESISTOGRAPH measurement to badly rotten trees. For those tree species that did not include living specimens, cross-sectional wood discs were collected from the stumps of felled healthy trees supplementally, and their ages were manually counted. Afterward, we estimated the regression equations for dominant tree species in the plots showing the relationship between the age after ingrowth (y) and the radius - 2.5 cm (w) of trees by polynomial equations of third order fitted to scatter diagrams. Table 2 shows the number of samples for each species and the estimated equations with coefficients of determination values. The appropriateness of the estimated equations was validated with RESISTOGRAPH data of the plots. These equations were applied to estimate the ages after ingrowth of badly rotten dead trees and living trees of all three plots that were impossible to estimate the ages from RESISTOGRAPH data and census data.

Table 2 Developed regression models based on RESISTOGRAPH measurements

Tree species	Number of tree samples	Equation	Coefficient of determination
Abies sachalinensis	51	y = −0.0000008w3 + 0.0005w2 + 0.1995w + 32.294	0.8379
Picea jezoensis	59	y = −0.000002w3 + 0.0008w2 + 0.2545w + 30.32	0.8506
Acer urkrunduense	41	y = −0.00001w3 + 0.0043w2 − 0.1416w + 41.991	0.7668
Tilia japonica	32	y = 0.000002w3 − 0.0012w2 + 0.5815w +1 3.189	0.7803

For the subsequent survival analysis, we classified the tree ages into the age classes of ten year length for each.

Survival Analysis of Target Trees

The survival data demonstrated that there were both dead and living trees in the study plots at the end of the observation periods. In these instances, the survival time is described as "censored". Right censoring happens when a target leaves the study before an event occurs, or the follow-up ends before an event occurs. On the other hand, "truncation" is used when we do not know whether events occurred before a specific time for targets.

The survival probability functions were developed based on the previous studies (Kalblfeisch and Prentice, 1980; Cox and Oakes, 1984; Crowder et al., 1991; Klein and Moeschberger, 1997; Kleinbaum and Klein, 2011; Fujikake, 2000; Hiroshima, 2014). The tree mortality is introduced as a stochastic variable (T) unique for each target tree. Here, T indicates the age class after ingrowth and is defined as a continuous variable that is based on the probability density function f (T). If a tree at the t-th age class dies, the conditional probability called the mortality rate (P_t) can be defined as follows:   

$\mathrm{Pr}(t-1<T\le t|T>t-1)=p_t$(1)

.

If a tree survives the t-th age class during the observation period, the conditional probability is defined as follows:   

$\text{Pr}(T>t|T>t-1)=1-p_t$(2)

.

By following the methods of Fujikake (2003), Hiroshima (2006), and Tiryana et al. (2011), and considering Eqs. 1 and 2, we can describe the likelihood function (L) of the observation as follows:   

$\begin{array}{c}L=\prod _t\mathrm{Pr}(t-1<T\le t|T>t-1{)}^{d_t}\mathrm{Pr}{(}^{T>t|}\\ =\prod _t{\left(\frac{{\int }_{t-1}^{t}f\left(T\right)dT}{{\int }_{t-1}^{\infty }f\left(T\right)dT}\right)}^{d_t}{\left(\frac{{\int }_t^{\infty }f\left(T\right)dT}{{\int }_{t-1}^{\infty }f\left(T\right)dT}\right)}^{a_t}\end{array}$(3)

where d_t is the number of dead trees and a_t is the number of surviving trees at the t-th age class during the observation period. By employing the probability density function, the mortality probability (q_t) of new ingrowth trees at the t-th age class is defined as follows:   

$\mathrm{Pr}(t-1<T\le t)=\underset{t-1}{\overset{t}{\int }}f\left(T\right)dT=q_t$(4)

.

In addition, the survival probability (r_t) at the t-th age class is defined as follows:   

$\mathrm{Pr}(T>t-1)=\underset{t-1}{\overset{\infty }{\int }}f\left(T\right)dT=1-\underset{1}{\overset{t-1}{\int }}f\left(T\right)dT=r_t$(5)

.

Therefore, based on Eq. 1, the mortality rate (P_t) can also be:   

$p_t=\frac{q_t}{r_t}$(6)

. Thus, Eq. 3 can be also expressed as:   

$\begin{array}{c}L=\prod _t{\left(\frac{{\int }_{t-1}^{t}f\left(T\right)dT}{{\int }_{t-1}^{\infty }f\left(T\right)dT}\right)}^{d_t}{\left(\frac{{\int }_t^{\infty }f\left(T\right)dT}{{\int }_{t-1}^{\infty }f\left(T\right)dT}\right)}^{a_t}\\ ={\prod _t\left(\frac{q_t}{r_t}\right)}^{d_t}\left(1-\frac{q_t}{r_t}\right)=\prod _t{p}_t{}^{d_t}{(}^1\end{array}$(7)

.

The maximum likelihood estimators of P_t can be calculated by solving the first-order derivation equation of Eq. 3 as shown in Eq. 8:   

$\stackrel{\wedge }{{\text{p}}_t}=\frac{d_t}{a_t+d_t}$(8)

.

Moreover, considering Eq. 5 and Eq. 6 lead to   

$r_t-r_{t+1}=q_t$(9)

, the survival probability can be converted into   

$r_{t+1}=r_t-q_t=r_t\left(1-\frac{q_t}{r_t}\right)=r_t\left(1-p_t\right)$(10)

.

This Eq. 10 can also be expressed as follows:   

$r_t=\prod _{k<t}\left(1-p_k\right)$(11)

.

Thus, the consistent estimator of r_t can be expressed as follows:   

$\stackrel{\wedge }{r_t}=\prod _{k<t}\left(1-\stackrel{\wedge }{p_k}\right)=\prod _{k<t}\left(1-\frac{d_k}{a_k+d_k}\right)$(12)

.

This consistent estimator is called the Kaplan–Meier estimates (Kaplan and Meier, 1958), which describe the distribution of survival probabilities. The Kaplan–Meier curves represent the non-parametric expressions of survival probability distributions throughout the study period.

A comparison of the survival probabilities can be made using the log-rank test (Mantel, 1966). This tests a null hypothesis (i.e., no significant difference in survival between consecutive periods in this study) and the expectation of an equal number of deaths (E) for each of the two groups. The observed (real) number of deaths is indicated with O in the following equation:   

$Logrankstatistic={\left(O-E\right)}^2/Var\left(O-E\right)$(13)

.

One of the advantages in survival analyses with covariates of age, i.e., age-based survival analyses, is to analyze and evaluate the maturity of cohorts considering age distributions, mean age, mean lifetime, etc. Particularly, the calculation of mean lifetime requires the above-estimated mortality probability like the following:   

$Meanlifetime=\sum t{q}_t/\sum q_t$(14)

.

In the case of this study, it is not proper to apply the non-parametric estimators of q_t in Eq. 14 when the dead trees are missing in some of the age classes because the denominator should be nearly equal to one. Otherwise, it is desirable to apply parametric survival analyses for mean lifetime calculations to smooth the stepwise non-parametric estimates. For that purpose, we apply the Weibull distribution for f(T) with the parameters of m and k:   

$f(t;m,k)=\frac{k}{m^k}{t}^{k-1}{e}^{-}{}^{{\left(\frac{t}{m}\right)}^k}$(15)

.

The mean and variance of Weibull distribution are $m\Gamma \left(1/k+1\right)$ and $m^2\left\{\Gamma \left(2/k+1\right)-{\Gamma }^2\left(1/k+1\right)\right\}$ respectively, and the former represents the mean lifetime. The parameters can also be estimated by the maximum likelihood method with L in Eq. 3.

RESULTS

Species Composition

The dead and living tree composition of three plots in three periods are listed in Table 3. In the plot 5203, the number of living trees for the period 1, 2 and 3 were 319, 295, and 260. Mortality percentage was 12.15%, 19.89%, and 18.51% for the three periods of 5203. Tilia japonica had the highest number of dead trees among all the species of 5203. Moreover, the highest percentage of 19.89% of tree mortalities among all three plots resulted from period 3 of 5203. For plot 5225, the number of living trees of three periods were 114, 193, and 189, respectively. The least tree mortality percentage of 5.80% among all the plots was from period 2 of 5225. In period 3 of 5225, Acer urkrunduense was the main dead species among all the species. Though Acer urkrunduense had a lower average tree height, it represented the highest number of tree mortalities during 5225. In the plot 5240, 207, 184, and 187 living trees were counted, and Abies sachalinensis and Picea jezoensis were the two major species. Dead tree percentages of 5240 were 6.91%, 9.39%, and 7.46% for the three periods accordingly. Considering differences in the species composition of living and dead trees among the three plots, we dealt with the three plots data separately in the following survival analyses.

Table 3 Living and dead species composition of each plot of all the periods (1989–2019)

Species	Plot 5203			Plot 5225			Plot 5240
	Total trees (dead trees)			Total trees (dead trees)			Total trees (dead trees)
	P1	P2	P3	P1	P2	P3	P1	P2	P3
Abies sachalinensis	91(7)	92(8)	94(20)	34(5)	45(2)	56(9)	50(5)	49(6)	63(4)
Picea jezoensis	62(10)	52(10)	42(10)	32(3)	31(1)	32(7)	42(0)	43(4)	40(3)
Picea glehnii	1(0)	1(1)	0(0)	1(0)	1(0)	1(0)	3(2)	1(1)	0(0)
Betura platyphylla	2(0)	3(0)	3(2)	6(1)	7(0)	9(2)	9(4)	5(0)	6(1)
Quercus crispula	0(0)	0(0)	0(0)	0(0)	0(0)	0(0)	2(0)	2(1)	1(0)
Ulmus laciniata	40(3)	38(7)	33(6)	7(1)	12(1)	11(1)	4(1)	3(0)	3(1)
Cercidiphyllumjaponicum	1(0)	1(0)	1(0)	0(0)	0(0)	0(0)	0(0)	0(0)	0(0)
Magnolia obovata	0(0)	0(0)	0(0)	1(1)	1(0)	1(1)	0(0)	0(0)	0(0)
Prunus maximowiczii	1(0)	1(0)	1(0)	0(0)	0(0)	0(0)	1(0)	1(0)	1(0)
prunus ssiori	4(0)	6(1)	6(0)	0(0)	1(0)	1(0)	7(0)	8(0)	8(1)
Sorbus commixta	3(0)	4(2)	2(0)	5(0)	8(2)	7(2)	29(3)	25(7)	23(6)
Sorbus alnifolia	1(0)	1(0)	1(0)	0(0)	0(0)	0(0)	2(0)	2(0)	2(0)
Phellodendron amurense	1(0)	1(0)	1(0)	0(0)	0(0)	0(0)	8(3)	5(1)	4(1)
Acer japonicum	10(2)	9(4)	11(2)	0(0)	2(0)	2(0)	9(2)	7(3)	4(0)
Acer ukurunduense	1(0)	1(1)	0(0)	28(8)	58(10)	69(19)	1(1)	0(0)	0(0)
Acer mono	6(0)	7(0)	8(1)	0(0)	0(0)	0(0)	9(0)	10(1)	9(2)
Acer mono var. myrii	19(2)	25(2)	27(1)	1(0)	1(0)	1(0)	12(1)	14(1)	13(1)
Tilia japonica	108(20)	110(35)	81(22)	21(3)	45(5)	45(8)	34(3)	32(7)	28(1)
Acanthopanax sciadophyloides	0(0)	0(0)	0(0)	0(0)	0(0)	0(0)	9(0)	9(1)	8(6)
Kalopanax pictus	8(0)	9(1)	8(0)	0(0)	0(0)	0(0)	1(0)	2(1)	1(0)
Aralia elata	0(0)	0(0)	1(1)	0(0)	0(0)	0(0)	0(0)	0(0)	0(0)
Fraxinus mandshuricavar. japonica	4(0)	6(0)	7(2)	2(2)	2(0)	2(0)	0(0)	0(0)	0(0)
Total trees (dead trees)	363(44)	367(72)	327(67)	138(24)	214(21)	237(48)	232(25)	218(34)	214(27)
Dead tree %	12.15	19.89	18.51	6.63	5.80	13.26	6.91	9.39	7.46

Age Class Distribution

Tree age class distributions are presented for three plots in three periods (Fig. 3). Note that the number of trees in the first age class consisted of the aggregation of new ingrowth trees found in census data in previous periods and that the number of trees in the t-th age classes were equal to the number of living trees in the (t-1)-th age classes in previous periods (t = 2, 3, 4,…). In plot 5203 (Fig. 3a), the highest number of new ingrowths can be seen from the 1st age class of period 1, while the lowest number of ingrowths from period 3. Most of the dead trees concentrated in the 1st, 2nd, and 3rd age classes through the periods (Fig. 3a). The distribution of 5203 showed many young trees and few old trees in the shape of the bimodal distribution in each period.

Fig. 3.

Age class distribution of living and dead trees of 3 plots over the periods.

In the plot 5225 (Fig. 3b), the highest number of new ingrowth trees were observed in the 1st age class of period 2. Period 3 had the highest number of dead trees among other periods. Periods 1 and 2 showed an exponential type of age class distribution while 2nd age class of period 3 had a remarkably higher number of trees as an exception. Thus, the exponential pattern of age class distribution presented in the plot 5225 except for the 2nd age class of period 3.

In the plot 5240 (Fig.3c), the 2nd and 7th age classes of period 1 had a higher number of trees among all the age classes of period 1. The 1st age class of period 2 had a lower number of new ingrowths compared with other periods. It affected the 2nd age class of period 3 as well. The highest number of dead trees came from the 3rd age class of period 2 (Fig. 3c). Therefore, all three periods almost showed the shape of exponential type though they had some exceptions.

Thus, the age class distribution patterns varied among the plots, so further survival analyses were carried out separately, i.e., we did not compare Kaplan–Meier curves among three plots in the following log-rank tests. Besides, due to the insufficient number of dead trees for each species, further analyses were carried out without considering the difference of tree species, i.e., considering only one group of trees.

Kaplan–Meier Curves

Figure 4 shows the estimated Kaplan–Meier curves in the three plots during three periods. In this study, we compared the consecutive periods of period 1 with 2 and period 2 with 3. The stepwise fluctuations were inherent to Kaplan–Meier curves. In all the three periods of the plot 5203 (Fig. 4a), the curves declined considerably before the 3rd age class as a result of many young tree mortalities and decreased stepwise in the age classes in which trees died. In Figure 4b of the plot 5225, compared with period 2, the curves declined considerably during periods 1 and 3 owing to the deaths of both young and mature trees. No decline was observed in the curve of period 2 for the 8th or higher age classes since no dead trees were observed after the 8th age class. The plot 5240 (Fig. 4c) showed a similar Kaplan–Meier curve distribution with the plot 5203 though the highest age classes were younger. Generally, in all plots and all periods, none of the curves reached zero under the existing oldest age class. Overall, the differences in mortality among periods of 5225 were statistically significant. The other two plots did not show a statically significant difference when considering consecutive periods as mentioned above (Table 4).

Fig. 4.

Kaplan–Meier curves showing the tree survival probability during the three study periods.

Table 4 Log-rank test results by comparing each period

Period	Plot 5203		Plot 5225		Plot 5240
	p-value	Chi-squared value	p-value	Chi-squared value	p-value	Chi-squared value
1989-1999 and 1999-2009	0.1833	1.7707	0.04088*	4.1857	0.2167	1.5258
1999-2009 and 2009-2019	0.5577	0.3436	0.0068*	7.3295	0.6289	0.2336

Significant relationships (p<0.05) are with *.

Parametric Estimates

Mean lifetime could not be calculated based on the estimated Kaplan–Meier estimates in the plots due to the lack of dead trees in several age classes (Fig. 3). Therefore, the Weibull distributions were applied to all three plots for smoothing the stepwise values of non-parametric estimates. Table 5 shows the estimated parameters of m and k and the mean and standard deviation of the Weibull distributions. Figure 5 illustrates the results of estimated parametric distributions of mortality probability, mortality rate, and survival probability in the period 3 of three plots and, for reference, all three periods of the plot 5225. The meanings of mortality rate, mortality probability, and survival probability were explained, for instances, by assuming 100 trees in period 3 of plot 5203 (Fig. 5a) like the followings: For the mortality probability of 0.05 at the 6th age class, if 100 new in-growth trees are produced now, 5 trees will die in the future at the 6th age class. For the mortality rates of 0.19 at the 6th age class, if 100 trees live in the 6th age class now, 19 trees will die at the 6th age class within this period. For the survival probability of 0.28 at the 6th age class, 28 trees will survive above the 6th age class among 100 new ingrowth trees.

Table 5 Probability distributions of the Weibull for each plot

	5203			5225			5240
Weibull parameter	Period 1	Period 2	Period 3	Period 1	Period 2	Period 3	Period 1	Period 2	Period 3
k	0.8767	0.8759	0.8566	1.0374	0.6042	0.9382	0.9825	0.9050	1.3522
m	8.3828	4.3157	3.8640	5.3051	20.7683	4.4234	8.6831	5.5414	7.8366
Mean	8.9517	4.6111	4.1836	5.2271	30.9680	4.5545	8.7497	5.8132	7.1840
Standard Deviation	10.24	5.28	4.90	5.04	53.97	4.86	8.91	6.43	5.37

m = scale parameter, k= shape parameter

Fig. 5.

Distribution of tree mortality probability, mortality rate, and survival probability based on Weibull distributions of three plots.

Figure 5 showed that the shapes of probability distributions were determined by the so-called shape parameter of k, while mean longevity was determined by both m and k of Table 5. In period 3, the plots 5203 and 5225 (Fig. 5a and 5c) showed decreasing trends of mortality probabilities and survival probabilities along with age classes. Among them, Figure 5a represented the highest morality probability and mortality rate values in the 1st age class and fell mostly in the 2nd age class and then decreased monotonously with an increase in age. Period 3 of plot 5225 (Fig. 5c) maintained a slightly constant mortality rate throughout all the age classes. Interestedly, 5240 (Fig. 5b) showed a different pattern for the mortality rate, which was an increasing trend along with the age class, which came from the different range of shape parameters k>1 for the period 3 of 5240 among all the plots. In addition, period 2 of 5225 (light gray dots in Fig. 5c) showed remarkably higher survival probabilities and the highest mean of 30.97 with the lowest k of 0.60 while all the other plots have the means within 4.18–8.95 (Table 5).

Calculation of Age-related Indicators

We calculated the mean lifetimes, mean stand ages, and mean biological lifetimes for the purpose of evaluating the matured states of study stands.

According to the Weibull estimates in Table 5, the mean value represents the mean lifetime after ingrowth, and the mean lifetimes of the plot 5203 were about 90, 46, and 42 years old for periods 1, 2, and 3, respectively (Table 6). The lowest mean lifetime among the three plots was 42 years from period 3 of this plot 5203. The mean lifetimes of the plot 5225 were 52, 310, and 46 years old for three periods. Consistent with the Kaplan–Meier curve in Figure 4b, the highest mean and standard deviation values were found in the period 2 of plot 5225. Similarly, the mean lifetimes of the plot 5240 were 87, 58, and 72 years old for three periods.

Table 6 Comparison of the mean lifetime, mean stand age, and mean biological lifetime of each plot

	Plot 5203			Plot 5225			Plot 5240
	Period 1	Period 2	Period 3	Period 1	Period 2	Period 3	Period 1	Period 2	Period 3
Mean lifetime	90	46	42	52	310	46	87	58	72
Mean stand age	33	39	46	45	31	35	43	49	51
Mean biological lifetime	163	162	161	169	167	166	169	168	165

The mean stand ages can be simply calculated as the weighted (i.e., the numbers of living trees in corresponding age classes are weights) means of tree ages in the study plots. These values ranged from 31 to 51 years old (Table 6).

The biological lifetimes of dominant species are reported as follows (Watanabe, 1970): Abies sachalinensis-127 years, Picea jezoensis-210 years, Acer spp.-178 years, and Tilia japonica-155 years. The mean biological lifetimes of study stands can be calculated as the weighted (i.e., the number of each species are weight) means of these biological lifetimes of the dominant species in the study plots. The mean biological lifespan of the study stands ranged from 161–169 years old (Table 6).

DISCUSSION

We discuss the matured state of stands from the viewpoints of age distribution, stable state of survival probability (KM curve), and comparison of age-related indicators referring to the previous studies.

Age Distribution

Matured forests can be viewed as aggregates of many patches in different stages of regeneration or cyclic succession, so that a wide variety of stand dynamics can lead to reverse-J shaped age distributions, particularly for natural stands composed of mixed tree species with different growth rates (Watt, 1947; Villalba and Veblen, 1997; Antos et al., 2000; Kuuluvainen et al., 2002; Worbes et al., 2003). In this study, the age distributions of all three plots declined gradually toward the older age classes resulting in two general shapes of the bimodal in the plot 5203 and the exponential (basically same meaning with the reverse J) in 5225 and 5240 except for several specific age classes. Hiroshima (2014) suggested that reverse J-shaped age structures would be maintained in enough matured stands assuming continuous new ingrowth and dead trees in the older age classes. Because the three study plots were classified as the stand type of "coniferous selective cutting with poor regeneration" in UTHF, some of the stands may not maintain current exponential distributions in the future owing to a lack of an adequate number of new ingrowths continuously. Considering these perspectives, moreover, Mosseler et al. (2003) discussed the critical aspects of old-growth forests and suggested that if new ingrowth trees were lacking, the stand moved to the old-growth stage with uneven multi-modal age structures as well as reverse J-shaped structures. Thus, the states of multi-modal age distributions could also be conceived as one characteristic of structural features in old-growth forests.

The Stable State of Survival Probability over Time

Hiroshima (2014) discussed the matured states of secondary natural forest stands and suggested that survival probability distributions, i.e., Kaplan–Meier curves, would converge into a stable state over time, referring to the generalized normal forest concept based on the Gentan probability theory (Suzuki, 1984) in enough matured stands such as natural forests in a climax stage of successions. Moreover, our findings suggested that the definition of "matured state" is based on two criteria; when there is no statistical difference among Kaplan–Meier curves over time and if the survival probabilities of old age classes reach nearly zero values. Therefore, the matured state of the forest stand can be evaluated based on the stable state of the survival probability over time.

In this study, significant differences were not found among Kaplan–Meier curves in consecutive two periods in the plot 5203 and 5240. Both plots showed no significant differences, but their KM curves shifted or fluctuated slightly from period to period. Particularly the time-series changes of 5203 were reasonable; with the increase of older dead trees in the later period (Fig. 3a), KM curves shifted downward in spite of no significant difference between each period (Fig. 4a), which looked certainly getting close to the stable state of a curve over time with the survival probability in the oldest age approaching to zero. Compared with the plot 5203, the curves of 5240 showed fluctuating in spite of no significant difference as well, which first shifted downward and then went upward. The KM curves of plot 5240, therefore, looked still on the way to converge (Fig. 4c). Incidentally, the KM curve of plot 5225 in period 2 showed the exceptional pattern (Fig. 4b) without older dead trees over 8th age classes (Fig. 3b).

Comparison of Mean Stand Age, Mean Biological Lifetime, and Mean Lifetime of Stands

Mosseler et al. (2003) summarized the aspects of old-growth forests like the followings: The mean age of predominant species reached almost half the maximum lifetime for the individual species, and some old trees were close to their maximum lifetime. In addition, when the stand remained in an old-growth stage without adequate new ingrowth continuously, the mean lifetime for the individual species would be getting higher, which was consistent with the suggestion by Hiroshima (2014). In our results, half of the mean biological lifetimes of study stands were around 80 years old and these 80 years old were still higher than the mean stand ages ranged 31 – 51 years old in all cases (Table 6). Meanwhile, some old trees were close to or even older than the mean biological stand lifetime of 160 years old in the study stands (Fig. 3), particularly the plot 5203 observed the maximum tree ages over 200 years old. The time-series changes in mean lifetimes in the study stands fluctuated, not showing the tendencies of getting older like Mosseler et al. (2003) insisted. Basically, the mean lifetime cannot be calculated correctly before the mortality probability converged in a stable manner; if there were not enough older dead trees in immature stands and survival probabilities remained higher values up to the maximum age classes, which simply led to the calculations resulted in a higher lifetime of stands. In this sense, if the plot 5203 was approaching to the stable state, the mean lifetime of around 40 years old in period 3 might be close to the intrinsic mean lifetime of the enough matured state of the stand.

Together with all the maturity level indicators, these facts could be summarized as follows: Firstly, the plot 5203 had the bimodal age class distribution, and 5225 and 5240 had the exponential shapes with a few exceptions of some age classes. These two kinds of age distributions were conceivable to be a matured state of stands to a certain extent and considering the inadequate new ingrowth continuously in these stands, the shape of multi-modal distribution would be more general in enough matured stands. In this sense, the plot 5203 might be reaching a more matured state than 5225 and 5240. Secondly, the Kaplan–Meier curves were not significantly different over the periods in the plots except for the period 2 of 5225. These two stands were conceivable to be a matured state of stands to a certain extent, and particularly, the KM curves of 5203 looked gradually approaching to the stable state of a curve over time compared with the 5240 showing fluctuated behavior of curves. In this sense, the plot 5203 might Thirdly, the mean stand ages were still lower than the half value of mean stand biological lifetime in all cases, while some old trees were close to or even older than the mean biological stand lifetime, particularly in 5203. The mean lifetime of stands fluctuated and did not show clear aging tendencies, which indicated the study stands did not yet reach the enough matured states though there was a possibility that the mean lifetime of plot 5203 in period 3 was relatively close to the intrinsic mean lifetime of the enough matured stand. In these senses, the plot 5203 might be reaching the more matured state than 5225 and 5240. To sum up all these three aspects, it was consistently implied that the study stands might not get enough matured yet, but the plot 5203 showed further progress toward the matured state than 5225 and 5240.

CONCLUSION

In this study, we applied age-based survival analyses to the natural forest stands and estimated both non-parametric and parametric estimates of tree mortalities. Moreover, we evaluate the maturity of stands by comparing multi-temporal metrics such as age distributions, survival probabilities, mean age, and mean lifetime of stands. The analyses implied that the study stands might not get enough matured yet, but some plots showed further progress toward the matured state than others.

The mortality patterns and survival probabilities reported in this study would constitute a valuable reference for future studies to understand the stand dynamics of natural forests coherent to the morality of individual tree populations. However, for a better understanding of tree mortality and for reducing uncertainty in estimated survival probabilities, further research is needed. For instance, survival probabilities considering species differences can be more useful to identify the over-matured trees to harvest in the SSMS of UTHF.

ACKNOWLEDGMENTS

We thank the technical staff of the UTHF for their logistic support throughout the fieldwork and their help with the tree census data.

LITERATURE CITED

• Alvi, M.A., McArt, D.G., Kelly, P., Fuchs, M.A., Alderdice, M., McCabe, C.M. and Loughrey, M.B. (2015) Comprehensive molecular pathology analysis of small bowel adenocarcinoma reveals novel targets with potential for clinical utility. Oncotarget 6 : 20863–20874
• Amateis, R.L., Burkhart, H.E. and Liu, J. (1997) Modeling survival in juvenile and mature loblolly pine plantations. For. Ecol. Manag. 90 : 51–58
• Ancarani, A., Di Mauro, C., Fratocchi, L., Orzes, G. and Sartor, M. (2015) Prior to reshoring: A duration analysis of foreign manufacturing ventures. Int. J. Prod. Econ. 169 : 141–155
• Antos, J.A., Parish, R. and Conley, K. (2000) Age structure and growth of the tree-seedling bank in subalpine spruce-fir forests of south-central British Columbia. Am. Midl. Nat. 143 : 342–354
• Buckley, N.E., Haddock, P., Simoes, R.D.M., Parkes, E., Irwin, G., Emmert-Streib, F. and Mullan, P. (2016) A BRCA1 deficient, NFκB driven immune signal predicts good outcome in triple negative breast cancer. Oncotarget 7 : 19884
• Cox, D.R. and Oakes, D. (1984) Analysis of survival data. (Vol. 21). Chapman and Hall, London, 212 pp
• Crowder, M.J., Kimber, A.C., Smith, R.L. and Sweeting, T.J. (1991) Statistical analysis of reliability data. Chapman and Hall, London: 1–11
• Eid, T. and Tuhus, E. (2001) Models for individual tree mortality in Norway. For. Ecol. Manag. 154 : 69–84
• Flewelling, J.W. and Monserud, R.A. (2002) Comparing methods for modelling tree mortality. In: Crookston, N.L. and Havis, R.N. (eds) Second Forest Vegetation Simulator Conference. USDA Forest Service, Rocky Mountain Research Station: 168-177
• Fujikake, I. (2000) An analyses of harvesting activities in a forest management entity using harvesting age distribution (tentative translation by author). Doctoral thesis, Kyoto University, Kyoto
• Fujikake, I. (2003) Estimation of Gentan probability based on the forest resource table. Proc. Inst. Stat. Math. 51 : 95–109
• Gross, S.R., O'Brien, B., Hu, C. and Kennedy, E.H. (2014) Rate of false conviction of criminal defendants who are sentenced to death. Proceedings of the National Academy of Sciences 111 : 7230–7235
• Hiroshima, T. (2006) Study on the methodology of calculating mean and variance of felling age in forest planning. Jpn. J. For. Plann. 40 : 139–149
• Hiroshima, T. (2014) Applying age-based mortality analysis to a natural forest stand in Japan. J. For. Res. 19 : 379–387
• Jung, E.Y., Baek, C. and Lee, J.D. (2012) Product survival analysis for the App Store. Marketing Letters 23 : 929–941
• Kalbfleisch, J.D. and Prentice, R.L. (1980) The statistical analysis of failure time data. NY John Wiley and Sons, New York, 360 pp
• Kaplan, E.L. and Meier, P. (1958) Nonparametric estimation from incomplete observations. J. Am. Stat. Assoc. 53 : 457–481
• Kato, R. (1952) The vegetation of the Tokyo University Forest in Hokkaido. Bull. Tokyo Univ. For. 43 : 1–18
• Klein, J.P. and Moeschberger, M.L. (1997) Survival analysis: Techniques for censored and truncated Data. Springer, New York: 1–87
• Kleinbaum, D.G. and Klein, M. (2011) Survival Analysis: A self-learning text. 3rd ed. Springer-Verlag, New York: 1–359
• Kuuluvainen, T., Mäki, J., Karjalainen, L. and Lehtonen, H. (2002) Tree age distributions in old-growth forest sites in Vienansalo wilderness, eastern Fennoscandia. Silva. Fenn. 36 : 169–184
• Mantel, N. (1966) Evaluation of survival data and two new rank order statistics arising in its consideration. Cancer Chemother. Rep. 50 : 163–170
• Monserud, R.A. (1976) Simulation of forest tree mortality. Forest Sci. 22 : 438–444
• Morse, B.W. and Kulman, H.M. (1984) Plantation white spruce mortality: estimates based on aerial photography and analysis using a life-table format. Can. J. For. Res. 14 : 195–200
• Mosseler, A., Lynds, J.A. and Major, J.E. (2003) Old-growth forests of the Acadian Forest Region. Environ. Rev. 11 : 47–77
• Nothdurft, A. (2013) Spatio-temporal prediction of tree mortality based on long-term sample plots, climate change scenarios and parametric frailty modeling. For. Ecol. Manage. 291 : 43–54
• Rinn, F., Schweingruber, F.H. and Schar, E. (1996) Resistograph and X-ray density charts of wood competitive evaluation if drill resistance profiles and X-ray density charts of different wood species. Holzforschung 50 : 303–311
• Suzuki, T. (1959) Yield prediction with the transition probability. Trans. Mtg. Kanto Branch Jpn. For. Soc. 11 : 36–38
• Suzuki, T. (1984) The Gentan probability, a model for the improvement of the normal forest concept and of forest planning. In: Nagumo, H., Konohira, Y., Kobayashi, S., Minowa, M., Nishikawa, K., Naito, K., Sweda, T., Amano, M. and Tanaka, K. (eds) Proceedings IUFRO symposium on forest management planning and managerial economics. The Organizing Committee of IUFRO Symposium on Forest Management Planning and Managerial Economics, Tokyo: 12–24
• Takahashi, N. (2001) The stand based forest management system: its theory and practice (rev.ed). Log Bee, Sapporo, 125 pp (in Japanese)
• Tatsuhara, S. and Minowa, M. (1988) Studies on the tansition of diameter distribution in coniferous plantations: An estimation of mortality probability by diameter class. Bull. Tokyo. Univ. For. 80 : 203–255
• Tiryana, T., Tatsuhara, S. and Shiraishi, N. (2011) Empirical models for estimating the stand biomass of teak plantations in Java, Indonesia. J. For. Plan. 16 : 35–44
• The University of Tokyo Hokkaido Forest (2019) http://www.uf.a.u-tokyo.ac.jp/files/gaiyo_hokkaido.pdf. (accessed 3 February 2020)
• Vanclay, J.K. (1995) Synthesis: growth models for tropical forests: a synthesis of models and methods. Forest Sci. 41 : 7–42
• Villalba, R. and Veblen, T.T. (1997) Regional patterns of tree population age structures in northern Patagonia: climatic and disturbance influences. J. Ecol. 85 : 113–124
• Watanabe, S. (1970) Hokkaido natural forestsucession pattern 1. Hokkaido For. 261 : 349–356
• Waters, W.E. (1969) The life table approach to analysis of insect impact. J. For. 6 : 300–304
• Watt, A.S. (1947) Pattern and process in the plant community. J. Ecol. 35 : 1–22
• Woodall, C.W., Grambsch, P.L. and Thomas, W. (2005) Applying survival analysis to a large-scale forest inventory for assessment of tree mortality in Minnesota. Ecol. Modell. 189 : 199–208
• Worbes, M., Staschel, R., Roloff, A. and Junk, W.J. (2003) Tree ring analysis reveals age structure, dynamics and wood production of a natural forest stand in Cameroon. For. Ecol. Manage. 173 : 105–123
• Wyckoff, P.H and Clark, J.S. (2000) Predicting tree mortality from diameter growth: a comparison of maximum likelihood and Bayesian approaches. Can. J. For. Res. 30 : 156–167