Age-Based Survival Analysis of Coniferous and Broad-Leaved Trees: A Case Study of Preserved Forests in Northern Japan

Abstract

Scientifically sound methods are essential to estimate the survival of trees, as they can substantially support sustainable management of natural forest resources. Tree mortality assessments have mainly been based on forest inventories and are mostly limited to planted forests; few studies have conducted age-based survival analyses in natural forests. We performed survival analyses of individual tree populations in natural forest stands to evaluate differences in the survival of two coniferous species (Abies sachalinensis (F. Schmidt) Mast. and Picea jezoensis var. microsperma) and all broad-leaved species. We used tree rings and census data from four preserved permanent plots in pan-mixed and sub-boreal natural forests obtained over 30 years (1989–2019). All living trees (diameter at breast height ≥ 5 cm in 1989) were targeted to identify tree ages using a Resistograph. Periodical tree age data, for a 10-year age class, were obtained during three consecutive observation periods. Mortality and recruitment changes were recorded to analyze multi-temporal age distributions and mean lifetimes. Non-parametric survival analyses revealed a multi-modal age distribution and exponential shapes. There were no significant differences among survival probabilities of species in different periods, except for broad-leaved species, which had longer mean lifetimes in each period than coniferous species. The estimated practical mean lifetime and diameter at breast height values of each coniferous and broad-leaved tree can be applied as an early identification system for trees likely to die to facilitate the Stand-based Silvicultural Management System of the University of Tokyo Hokkaido Forest. However, the survival probabilities estimated in this study should be used carefully in long-term forest dynamic predictions because the analysis did not include the effects of catastrophic disturbances, which might significantly influence forests. The mortality patterns and survival probabilities reported in this study are valuable for understanding the stand dynamics of natural forests associated with the mortality of individual tree populations.

1. Introduction

Changes in the survival probability of forests have severe environmental [1] and economic consequences [2] and can influence forest management decisions [3,4]. Hence, long-term tree mortality studies are needed [5,6,7,8] considering the threat of global climate change, which is forecasted to increase tree mortality in forests worldwide [9,10,11,12,13]. Because of the difficulty in identifying the exact ages of trees, researchers mainly use diameter at breast height (DBH), dominant height, basal area, growth rate, or competition index, avoiding age-based measures to determine mortality probabilities [14] Popular methods for identifying tree mortality are logistic models [15], Weibull [16], the Gamma [17], Richard’s function [18], and exponential or probit models [19]. Tree age can be used for the accurate prediction of tree mortality [20,21,22]. Generally, the lifespan of a tree is unknown; unlike logistic regressions, the right censoring and left truncation approaches of survival analysis can handle tree age data [14,23]. Several survival analysis techniques have been suggested to determine forest mortality [24]; however, these are limited to even-aged forest stands [25,26,27]. Woodall et al. [14] carried out a survival analysis based on DBH. More recently, survival analysis of natural forests has been successfully applied based on inventory measurements (i.e., [28,29,30]). Nothdurft [28] provided spatiotemporal predictions of tree mortality based on parametric frailty modeling. Moreover, survival models have been fitted based on inventory data by omitting the first 20 years of a tree’s life [31]. However, these studies emphasized the importance of determining exact tree ages for the survival analysis of forests [14,29,30].

The principles underlying the approach presented in this paper are related to the idea presented by Hiroshima [32], which involved more flexible techniques to find tree age for survival analysis based on non-parametric and parametric approaches. The presented method is based on annual tree ring measurements obtained using a semi-nondestructive device, as well as periodic inventory data of secondary natural forests in Japan. However, previous studies focusing on the age-based survival analysis based on uneven-aged natural forests are comparatively rare, because they require unique methods to detect tree age and conduct long-term, large-scale observations [33,34].

In this study, we used a device to determine tree age at the long-term research site in the preserved area in the University of Tokyo Hokkaido forest (UTHF), which is considered ideal for performing investigations similar to the present one, because plots in this area do not undergo any management following natural disturbance events [35]. Overall, the purpose of this study was to perform a survival analysis of individual tree populations in natural forest stands of the UTHF, Japan, and estimate and compare the age distribution of two major coniferous species (Abies sachalinensis (F. Schmidt) Mast. and Picea jezoensis var. microsperma) and broad-leaved species. We applied survival analysis techniques to each of three observational periods in a 10-year age class to determine the non-parametric survival probabilities of the species. Specifically, we revealed the multi-temporal age distribution changes and mean lifetime changes in two groups of trees (coniferous and broad-leaves) during the three observational periods. Survival probability changes over the periods can be applied to determine the tree harvesting strategy of the UTHF.

2. Materials and Methods

2.1. Inventory and Site Data

Retrospective inventory and site data were derived from UTHF, located in Furano (central Hokkaido, Japan). This site experiences a mean annual temperature and precipitation of 6.4 °C and 1297 mm, respectively [35]. Dominant stratified soils include brown forest soil, dark brown forest soil, black soil, and podzol [36]. Mixed conifer–hardwood forest, typical of the cool temperate zone, covers most of the site area. Tree species commonly found in this forest area are Fraxinus mandshurica, Ulmus davidiana var. japonica, Alnus hirsuta, and Salix spp. in deciduous swamp forests at lower elevations of less than 300 m; a coniferous and broad-leaved mixed forest dominated by A. sachalinensis at middle elevations (300–600 m), scattered forests mixed with P. jezoensis, Picea glehnii, and Betula ermanii at upper elevations (800–1200 m), and alpine vegetation (e.g., Pinus pumila) in the upper forest limit (>1200 m).

The Stand-based Silvicultural Management System (SSMS), in other words, a natural forest management system based on selected cutting and natural regeneration, is currently being employed extensively in UTHF except for research plots [37,38]. In SSMS, 10–17% of the stand volume is harvested by single-tree selection, with a cutting cycle of 15–20 years, by removing defective (e.g., diseased, senescent, non-vigorous, and twisted) and over-matured trees. This system helps maintain tree health and productivity in the stand and controls the stand composition [37]. To determine cutting rates, permanent plots and long-term ecological research plots have been established and periodically assessed. Among these permanent plots, 25 are located in the preserved area, ranging in size from 0.04 to 2.25 ha, and have an elevation between 380 and 1290 m. Within these plots, DBH measurements of all trees with a DBH ≥ 5 cm are regularly performed by UTHF staff; in most cases, 5-year interval and census data are available for the last five decades. Other than these periodical measurements, no human intervention has occurred in this preserved area for several decades.

2.2. Survival Data

Survival data were derived from periodic surveys on preserved permanent sample plots. We selected four plots, ranging in size from 0.04 to 2.25 ha, located at an elevation range of 570 to 690 m with similar slope aspects and slope angles (Figure 1;

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 (°)
5203	0.40	580	327	Southwest	15
5224	0.25	570	150	Southwest	18
5225	0.25	690	237	Southwest	20
5240	0.25	600	214	Southwest	25

). The main soil types found in the plots are brown forest and podzolic soils [35]. In addition, the stands in the four plots are all classified as “coniferous selective cutting with poor regeneration” in UTHF, where no continuous sufficient and new in-growth trees are expected. All four plots were within close proximity and had similar species composition. The typical vegetation of the plots was coniferous and broad-leaved mixed forests dominated by A. sachalinensis, P. jezoensis, Acer spp., and Tilia spp. Therefore, the four plots were aggregated in further analyses on tree survival.

We used tree census data (species, DBH, living or dead status, cause of death, etc.) of the plots, obtained between 1989 and 2019. Within this period, we set three observation periods of 1989–1999 (period 1), 1999–2009 (period 2), and 2009–2019 (period 3). We detected the number of tree rings at breast height (1.3 m) in 2019 using a semi-nondestructive device, RESISTOGRAPH^® [39]. All target trees were living and had a DBH ≥ 5 cm in 1989; some of these trees died by 2019. Live/dead trees and new in-growth trees were counted at the end of each observation period. The connected RESISTOGRAPH (Heidelberg, Germany) data logger recorded measurements, which were transferred to a computer for further analyses. The field measurement tree ring data (radius at breast height −2.5 cm) in 2019 were extracted via the DECOM^TM software (Philadelphia, PA, USA), which is used for annual tree ring detection. The “radius at breast height −2.5 cm” reveals the tree age after in-growth if we set DBH = 5 cm as an in-growth border.

Figure 2 shows the annual rings in one living A. sachalinensis tree detected using the DECOM^TM software. RESISTOGRAPH measurements were ineffective for severely rotten and center-decayed trees, for which we established the regression equation between the age after in-growth (y) and the radius −2.5 cm (w) [40]. Next, the age after in-growth was estimated by inputting the radius data into simple three-dimensional equations fitted to a scatter diagram (

Table 2. Regression models developed based on RESISTOGRAPH measurements.

Tree Species	No. of Tree Samples	Equation	Coefficient of Determination
A. sachalinensis	51	y = $-$0.0000008w3 + 0.0005w2 + 0.1995w + 32.294	0.8379
P. jezoensis	59	y = $-$0.000002w3 + 0.0008w2 + 0.2545w + 30.32	0.8506
A. ukurunduense	41	y = $-$0.00001w3 + 0.0043w2 + 0.1416w + 41.991	0.7668
T. japonica	32	y = 0.000002w3 − 0.0012w2 + 0.5815w + 13.189	0.7803
Other broad-leaved species	45	y = $-$0.000005w3 + 0.0023 w2 + 0.0088 w + 42.386	0.8592

). This method was followed for 20% of the sample trees.

Table 3. Methodologies used to determine the age of each tree species.

Species	Census Data		Direct AResistograph Measurements		Developed ARegression Models		Total No. of Trees
	Living Trees	Dead Trees	Living Trees	Dead Trees	Living Trees	Dead Trees	Living Trees	Dead Trees
A. sachalinensis	113	9	58	16	48	16	219	41
P. jezoensis	18	4	55	11	49	6	117	21
Broad-leaved trees	289	67	76	21	55	10	423	98
Percentage	45.70	8.71	20.57	5.22	16.54	3.48

further illustrates the method applied to determine the in-growth years of trees.

Survival models were built for two major coniferous species present in UTHF, A. sachalinensis and P. jezoensis. All broad-leaved trees were considered collectively to develop survival models. The model data of period 3 were for 260 A. sachalinensis, 138 P. jezoensis, and 521 broad-leaved trees (see Table 3 for the number of dead and living trees in period 3 included in this study). The ages of target trees in period 3 were first identified using the above-mentioned methods. The ages in periods 1 and 2 were calculated by deducting 10 and 20 years from those in period 3. Finally, tree ages were classified, with 10 years per age class.

2.3. Fundamentals of Survival Analysis of Major Species

Survival probability functions were the essential elements of this analysis, as they reflect the overall performance of the major species; they were developed based on previous studies [20,22,32,41,42,43]. The survival data demonstrated that both dead and living trees were present in the study plots at the end of the observation period. When we carry out survival analysis, the issue of missing data appears most of the time. Therefore, the censoring and truncation techniques must be incorporated to overcome this issue. Right censoring occurs when a target leaves the study before an event occurs or the follow-up ends before mortality happens. In contrast, left truncation occurs when an object is not observed from the beginning of the study, but rather enters the study at a later time in the observation period [20,44]. Furthermore, if target trees survived at the t-th age class in the observation period, then observed mortalities were truncated and censored at both beginning and end of the observation period.

Mortality rates ($p_t$) were calculated using age class after in-growth (T), according to Equation (1):

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

(1)

where $p_t$ is the mortality rate in the t-th age class.

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

$$\mathrm{Pr}(\left.T>t\right|T>t-1)=1-p_t.$$

(2)

Following previous methods [45,46,47] and considering Equations (1) and (2), we described the likelihood function (L) of the observation as follows:

$$L={\displaystyle \prod _t\mathrm{Pr}(}\left.t-1<T\le t\right|T>t-1{)}^{d_t}\mathrm{Pr}{(\left.T>t\right|T>t-1)}^{a_t}={{\displaystyle \prod _t{p}_t^{d_t}\left(1-p_t\right)}}^{a_t},$$

(3)

where d_t is the number of dead trees and a_t is the number of surviving trees in the t-th age class during the observation period.

The mortality probability ($q_t$) of a new in-growth tree in the t-th age class was defined as follows:

$$\mathrm{Pr}(t-1<T\le t)=q_t.$$

(4)

The survival probability ($r_t$) in the t-th age class was defined as follows:

$$\mathrm{Pr}(T>t-1)=r_t.$$

(5)

Therefore, based on Equation (1), the mortality rate ($p_t$) can be expressed as follows:

$$p_t=\frac{q_t}{r_t}.$$

(6)

The maximum likelihood estimators of $p_t$ can be calculated by the first-order derivation of Equation (3), as shown in Equation (7):

$$\stackrel{\wedge }{{\mathrm{p}}_t}=\frac{d_t}{a_t+d_t}.$$

(7)

Considering that Equations (5) and (6) can be combined as follows:

$$r_t-r_{t+1}=q_t,$$

(8)

the survival function 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).$$

(9)

Equation (9) can also be expressed as follows:

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

(10)

The survival function is described by the following formula:

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

(11)

where $\stackrel{\wedge }{r_t}$ represents the survival function and $d_k/(a_k+d_k)$ represents the hazard function in the t-th age class.

This $\stackrel{\wedge }{r_t}$ is the Kaplan–Meier estimate, which is an important tool for analyzing censored data [23]. Survival analysis is performed to describe the distribution of tree mortality using Kaplan–Meier estimates in the form of step-wise curves [23].

It is meaningful to assess whether there are differences in survival (or cumulative incidence of the event) among different groups. For this purpose, there are several tests available to compare survival among independent groups.

Survival probabilities can be compared using the log-rank test [48], which 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) in each of the two groups of each species. The observed (i.e., real) number of deaths is indicated by O in the following equation:

$$\mathrm{Log}\mathrm{rank}\mathrm{statistic}={\left(O-E\right)}^2/\mathrm{Var}\left(O-E\right).$$

(12)

Wilcoxon test [49] can be used as an alternative to the log-rank test. It emphasizes the information at the beginning of the survival curve where the number at risk is significant, allowing early failures to receive more weight than later failures:

$$\mathrm{Wilcoxon}\mathrm{test}\mathrm{statistics}=\frac{{\left({\displaystyle \sum _{f}w(t_{(f)})(O-E)}\right)}^2}{\left({\displaystyle \sum _{j}w(t_{(f)})(O-E)}\right)}.$$

(13)

It weights the observed minus expected score at the time t_f by the number at risk, w(t_f) overall groups at time t_f.

Mean lifetime, which is also expressed as mean longevity or mean lifespan, can be considered as the area under survival curve. It is mathematically calculated as the weighted sum of the age and estimated mortality probability described previously herein and can be expressed as follows:

$$\mathrm{Mean}\mathrm{lifetime}={\displaystyle \sum t{q}_t}/{\displaystyle \sum q_t}.$$

(14)

It is meaningful to calculate not only mean lifetime but also “practical” mean lifetime from the forest management perspective based on the survival estimates for the two groups of trees (coniferous and broad-leaved). Common mean lifetime was derived by considering all the age classes whereas practical mean lifetime used only specific age class following the Kaplan–Meier estimates. The particular point of the practical mean lifetime is that it has avoided the drastic decrease of survival probabilities in younger age classes, which commonly happen in natural forest stands due to suppression. We intended to explore the possibility of using the practical mean lifetime for kinds of harvesting standards in tree marking process, so that too short mean lifetime was inconvenient. Finally, these practical mean lifetime values were converted into DBH values by using the equations of Table 2 to facilitate practical application in the tree marking process of SSMS.

3. Results

3.1. Temporal Forest Structure Changes

The dominant species with relatively large sample sizes are listed in

Table 4. Summary of the dominant tree species in the four plots in the three observation periods.

Species	No. of Trees (%)
	Period 1		Period 2		Period 3
	Living Trees	Dead Trees	Living Trees	Dead Trees	Living Trees	Dead Trees
Conifer
A. sachalinensis	201 (21.90)	22 (2.40)	215 (22.16)	22 (2.27)	219 (23.60)	41 (4.42)
P. jezoensis	148 (16.12)	14 (1.53)	135 (13.92)	16 (1.65)	117 (12.61)	22 (2.37)
Other conifers	4 (0.44)	2 (0.22)	1 (0.10)	3 (0.31)	1 (0.11)	0 (0.00)
Total conifer	353 (38.45)	38 (4.14)	351 (36.19)	41 (4.23)	337 (36.31)	63 (6.79)
Broad-leaved
Ulmus laciniata	56 (6.10)	8 (0.87)	52 (5.36)	10 (1.03)	46 (4.96)	8 (0.86)
Sorbus commixta	35 (3.81)	3 (0.33)	27 (2.78)	11 (1.13)	25 (2.69)	8 (0.86)
Acer ukurunduense	21 (2.29)	9 (0.98)	48 (4.95)	11 (1.13)	50 (5.39)	19 (2.05)
Acer mono var. myrii	56 (6.10)	2 (0.22)	57 (5.88)	6 (0.62)	54 (5.82)	7 (0.75)
Tilia japonica	174 (18.95)	29 (3.16)	171 (17.63)	52 (5.36)	153 (16.49)	32 (3.45)
Other broad-leaved	115 (12.53)	19 (2.07)	113 (11.65)	20 (2.06)	101 (10.88)	25 (2.69)
Total broad-leaved	457 (49.78)	70 (7.63)	468 (48.25)	110 (11.34)	429 (46.23)	99 (10.67)
Total	810 (88.24)	108 (11.76)	819 (84.43)	151 (15.57)	766 (82.54)	162 (17.46)

No. of trees (%) = Number of trees in each period of major species and percentage of it within brackets.

. The demographic parameters of all the species were analyzed; however, only the major species in each period are presented in Table 4. In 1989 (period 1), there were 28 species, and the majority of living trees were broad-leaved species (49.78%). The percentage of coniferous species was highest in period 1, with the highest number of trees being that of A. sachalinensis. Following the same tendency of period 1, periods 2 and 3 also had higher percentages of broad-leaved species than of coniferous species. In 1999 (period 2), 27 species were sampled in the four permanent sample plots, and the dominant living species were A. sachalinensis (22.16%) and Tilia japonica (17.63%). In period 3, the dominant species among living trees was A. sachalinensis (23.60%). The number of dead coniferous trees increased over the three periods, and most of them were A. sachalinensis. Tilia japonica was the most common species among dead broad-leaved trees, and it showed an increasing trend, except in period 3.

A. sachalinensis and P. jezoensis accounted for approximately 99% of the coniferous species in each period, and the results for these species, with relatively large sample sizes, were considered for further analysis. Due to the insufficient number of dead trees of each broad-leaved tree species, all broad-leaved trees were combined in further analyses.

3.2. Changes in Age Distributions by Species Groups

In the three periods (periods 1–3), tree age-class distributions, including dead and living trees by three species groups, are presented in Figure 3. New in-growths found in the census data were aggregated in the 1st age class of each period.

In Figure 3a, the highest number of new in-growths can be seen in the 1st age class of period 3, whereas the lowest number of in-growths can be seen in period 1. The dead A. sachalinensis trees were distributed almost equally across all age classes. The distribution of A. sachalinensis trees showed an exponential shape, with many young trees and few old trees, in each period, with a few exceptions in period 1; for example, the 3rd age class of period 1 had a higher number of trees than the 2nd age class.

For P. jezoensis in Figure 3b, the largest number of new in-growth trees was observed in the 1st age class of period 1 among the three periods. Dead trees were distributed equally across those present in most age classes. P. jezoensis had a relatively low number of new in-growths in each period compared with A. sachalinensis, which might lead to a multi-modal age distribution, different from that of A. sachalinensis.

Figure 3c shows the distribution of all broad-leaved trees found in each period, as mentioned above. The largest number of new in-growths can be seen in the 1st age class of period 1. The 1st age class of period 3 had a lower number of new in-growths than those in other periods. Dead trees were distributed among almost all the age classes across the periods. Periods 1 and 2 exhibited an exponential age class distribution; however, the 1st age class of period 3 had a relatively low number of new in-growth trees. Thus, an exponential age class distribution was observed in broad-leaved species, except in the 1st age class of period 3.

3.3. Kaplan–Meier Curve Differences among the Periods

Figure 4 shows the estimated Kaplan–Meier curves for the three major species over the three observation periods. Consecutive periods (periods 1–3) are compared here. For A. sachalinensis, Figure 4a shows similar curves over the age classes in periods 1 and 2, with constant survival probabilities after age classes 15 and 14. In contrast, period 3 had a decreasing survival probability throughout age classes 1 to 20. The Kaplan–Meier curves for P. jezoensis show three different patterns for the three periods in Figure 4b. Remarkably, period 1 had higher survival probabilities throughout all age classes. The curve of period 3 shows a lower survival probability than that of period 2 after the 10th age class. The curves of broad-leaved species in Figure 4c show similar distributions, though the values remained slightly higher for period 1 than for the others after the 10th age class. In the results of the Log-rank and Wilcoxon statistical tests (

Table 5. Log-rank test results comparing the three observational periods.

Period	Tree Species
	A. sachalinensis				P. jezoensis				All Brzoad-Leaved Species
	Log-Rank Test		Wilcoxon Test		Log-Rank Test		Wilcoxon Test		Log-Rank Test		Wilcoxon Test
	p-Value	Chi-Square Value	p-Value	Chi-Square Value	p-Value	Chi-Square Value	p-Value	Chi-Square Value	p-Value	Chi-Square Value	p-Value	Chi-Square Value
Periods 1 and 2	0.9161	0.0111	0.6542	0.2007	0.5023	0.4502	0.3674	0.8124	0.0194 *	5.4637	0.0772	3.1236
Periods 2 and 3	0.1034	2.6524	0.5685	0.3252	0.1538	2.0337	0.2549	1.2961	0.2691	1.2215	0.0610	3.5101

* Indicates a significant difference (p < 0.05).

), the differences in mortality of broad-leaved species between periods 1 and 2 were only statistically significant based on the log-rank test (p-value = 0.0194) Coniferous species did not show statistically significant differences, in spite of the differences in their appearance, such as in the case of P. jezoensis.

3.4. Kaplan–Meier Curve Differences by Species

Figure 5 shows the Kaplan–Meier curves of the two major conifer species, A. sachalinensis and P. jezoensis, for the three periods. In addition,

Table 6. Log-rank test results of Abies sachalinensis and Picea jezoensis in each period.

	A. sachalinensis and P. jezoensis
	Period 1		Period 2		Period 3
Statistical Test	p-Value	Chi-Square Value	p-Value	Chi-Square Value	p-Value	Chi-Square Value
Log-rank test	0.4985	0.4581	0.5751	0.3142	0.9959	0.0000
Wilcoxon test	0.4827	0.4928	0.5720	0.3193	0.8692	0.0271

shows the results of statistical tests for each period. There was no statistically significant difference for any combination. Therefore, we combined the two conifer species to determine species group differences in the next step.

Figure 6 shows the estimated Kaplan–Meier curves for the new two groups of all conifers and all broad-leaved species over the three observation periods. For period 1, Figure 6a shows that the curve of broad-leaved trees declined considerably more than that of the conifers due to the high mortalities of young and old trees. No decline was observed in the curve of coniferous trees for age class 15 or higher because no dead trees were observed after these age classes. For period 2, Figure 6b shows different Kaplan–Meier curves from those of period 1; broad-leaved trees reached nearly zero at older age classes, whereas conifers remained constant after age class 14, as there were no dead trees after that point. For period 3, Figure 6c shows that both of the curves declined constantly because of young and older tree deaths across age classes. Notably, both curves reached almost zero for the oldest existing age class.

Table 7. Log-rank test and Wilcoxon test comparing the conifer and broad-leaved trees in each period.

	Conifer and broad-Leaved Species
	Period 1		Period 2		Period 3
Statistical Tzst	p-Value	Chi-Square Value	p-Value	Chi-Square Value	p-Value	Chi-Square Value
Log-rank test	0.1614	1.9607	0.0097 *	6.6870	0.0381 *	4.3018
Wilcoxon test	0.0386 *	4.2797	0.0463 *	3.9719	0.3035	1.0586

* indicates a significant difference (p < 0.05).

shows the results of statistical tests. Overall, differences in mortality between conifer and broad-leaved trees in each period were statistically significant in either log-rank or Wilcoxon tests.

3.5. Mean Lifetime over Three Periods

The (common) mean lifetimes after in-growth were calculated based on non-parametric estimates to investigate the temporal changes between coniferous and broad-leaved species (

Table 8. Mean lifetime and standard deviation (in parenthesis) of two groups of trees (coniferous and broad-leaves) in each period.

Period	Conifer Species	Broad-Leaved Species
Period 1	39 (4)	44 (5)
Period 2	36 (4)	36 (5)
Period 3	43 (5)	42 (5)

). The mean lifetimes of all conifer species were 39, 36, and 43 years for periods 1–3, respectively, and those of all broad-leaved species were 44, 36, and 42 years, respectively. Table 7 shows the results of statistical tests. Overall, differences in mortality between conifers and broad-leaved trees in each period were statistically significant in either log-rank or Wilcoxon tests.

Next, the “practical” mean lifetimes were calculated for period 3 only. We calculated this for period 3 only, because the Kaplan–Meier curves for that period appeared to almost converge to time-stable states, with probabilities reaching almost zero (Figure 6) and no significant difference between periods 2 and 3 in either tree category (Table 5). In addition, because there was no significant difference between the two major coniferous species (Table 6), it was suitable to derive one practical mean lifetime value for all coniferous species. Thus, the practical mean lifetime of all conifer groups was calculated as 72 years when avoiding the first five age classes, illustrated as line (ii) in Figure 6c (

Table 9. The predicted practical mean lifetime and corresponding diameter at breast height value of period 3.

Species	Practical Mean Lifetime in Years (SD)	Practical Mean DBH in cm (SD)
Conifer species	72 (7)	31 (7)
Broad-leaved species	80 (5)	36 (8)

Diameter at breast = DBH; standard deviation = SD.

). For all broad-leaved species, it was 80 years when avoiding trees of the first three age classes illustrated as line (i) in Figure 6c (Table 9). These two lines were identified as the age classes up to which tree mortality due to suppression occurred, and survival probability curves drastically declined in a single step-wise manner.

Finally, the practical mean lifetimes were converted to the corresponding DBH of 31 cm and 36 cm using the DBH-age equations developed (Table 9).

4. Discussion

4.1. Temporal Changes in Age Distribution and Species Composition

In Figure 3, the age distributions of A. sachalinensis and broad-leaved species show exponential shapes, with few outliers in some periods, whereas P. jezoensis shows an irregular multi-modal shape in each period. Remarkably, fewer new in-growths of coniferous species were found compared with those of broad-leaved species. Considering these matters, it appears that A. sachalinensis might not keep the current exponential shape and change to the multi-modal shape in the upcoming periods. In addition, larger numbers of dead trees and new in-growth trees were observed in broad-leaved species compared with conifer species over these periods. These results could not fully support the suggestion that the dominance of coniferous species would decrease further with climate change in northern Hokkaido, as found by Hiura et al. [50].

4.2. Survival Analyses over the Observational Periods

Regarding period differences, no significant differences were found between the Kaplan–Meier curves of two consecutive periods for the coniferous species, although the Kaplan–Meier curves declined along with the periods. Broad-leaved species in period 1 had a significantly greater survival probability than those in period 2. The only significant difference in the Kaplan–Meier curves was found between periods 1 and 2 for broad-leaved species, which might have been caused by the decrease in mortality in period 1. These survival probability trends correspond to the findings of Sato [51] and Hiura et al. (1996) [52]. Moreover, Kaplan–Meier curves between periods 2 and 3 were not significantly different for either tree category. The findings of Wijenayake and Hiroshima [40] suggested that the stand could be considered as an almost matured state when there was no statistical difference among Kaplan–Meier curves over time and the survival probabilities of old age classes reached nearly zero. Along with this, the survival probability values of period 3 can be applied to predict future age distributions.

Regarding species differences, the survival of coniferous and broad-leaved trees was statistically different over the three periods. Nakashizuka [53] also reported that the growth and survival of coniferous and broad-leaved species reflected large variabilities along with time. Therefore, it is reasonable to differentiate tree groups as coniferous and broad-leaved trees in further analyses.

4.3. Mean Lifetime Values and Harvesting Decisions

Selvin [54] discussed the importance of analyzing age at the time of death to provide more accurate lifetime estimations on survival probability. Yamamoto [55] indicated that natural mortalities can be minimized based on SSMS, although SSMS could fail if it does not adequately identify over-matured trees prior to decay. The early identification of over-matured trees may improve the SSMS practices, as well as sustainable natural forest management. Based on this belief, selection cutting can be planned and implemented to improve the quality and quantity of crop trees. The early identification of over-matured trees facilitates the maintenance of stands in continuously healthy and productive states [37].

The mean lifetimes did not differ greatly between conifers and broad-leaves; both were approximately only 40 years (Table 8). Instead of these mean lifetimes, the developed practical mean lifetimes might be applied for harvesting decisions in SSMS because they exclude juvenile tree mortalities due to suppression. These values were estimated based on the step-down width of each Kaplan–Meier curve and were inherent to species. In natural forest management, it is much more sensible in terms of DBH values rather than age values to facilitate harvesting decisions. The estimated DBH values of coniferous and broad-leaved species based on practical mean lifetime were 31 cm and 36 cm, respectively. These DBH values could be supplementary tools for the prediction of over-mature trees in SSMS. For example, with broad-leaved trees, it is best to pay attention when their DBH reaches 36 cm. This attention allows for selecting the harvestable trees in each cutting period by maintaining vigorous trees in the stand. Therefore, these developed practical mean lifetime and DBH values would be a supplementary tool in the practice of SSMS for early identification of the trees likely to die in the near future.

Based on the principles of SSMS, forest stands must be maintained in the pre-climax stage, which enables a high growth rate at the stand level [56]. Forest managers need a comprehensive understanding of natural stand development processes when designing silvicultural systems that integrate ecological and economic objectives [57]. Foresters usually manage stands that are developing as part of secondary succession. The significant tree species in terms of economic importance are commonly part of numerous seral stages towards the climax. Therefore, foresters manage their forests by controlling the tendency of that population to move toward a climax species forest [58]. By following the DBH, which is based on practical mean lifetime, stands can be maintained in the pre-climax successional stage.

The present analysis might not precisely reflect the mean lifetime in the future because the analysis primarily relied on non-parametric analyses. Parametric analysis approaches such as Weibull and Gamma could facilitate more precise future predictions, considering species differences, to identify likely dead trees for harvesting in SSMS. Therefore, in future research models, parametric analyses could be applied to predict age distributions.

5. Conclusions

This study showed the age class distribution of two coniferous tree species and all broad-leaved species at the long-term research site in the preserved area in UTHF over three observation periods and estimated survival probabilities using non-parametric methods. The resulted mean lifetime values were included suppressed young tree mortalities as well. Therefore, “practical” mean lifetimes were introduced by avoiding the drastic decrease in survival probabilities of Kaplan–Meier estimates. Furthermore, the developed values were converted to DBH values to enhance the practical use in the field. For example, when preparing the management plan, foresters should pay special care for the broad-leaved trees where DBH is 36 cm or more, and 31 cm for coniferous trees. Therefore, these developed practical mean lifetime and DBH values would be a supplementary tool in the practice of SSMS for early identification of the trees likely to die in the near future. Moreover, 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 associated with the mortality of individual tree populations.