Prediction of Tree Age Distribution Based on Survival Analysis in Natural Forests: A Case Study of Preserved Permanent Plots in the University of Tokyo Hokkaido Forest, Northern Japan

: In forests, tree mortality is strongly determined by complex interactions between multiple biotic and abiotic factors, and analysis of tree mortality is widely implemented in forest management. However, age-based tree mortality remains poorly evaluated quantitatively at the stand scale of the uneven-aged forests. The objective of this study is to predict the age distribution of living and dead trees based on survival analyses. We used a combination of tree-ring and census data from the two preserved permanent plots in the University of Tokyo Hokkaido Forest in pan-mixed and sub-boreal natural forests, Hokkaido, northern Japan, to derive site-specific survival models. All the living trees (diameter at breast height ≥ 5 cm in 2009) were targeted to identify tree ages using a RESISTOGRAPH, a semi non-destructive device. Periodical tree age data with a 10-year age class were used during the observation periods of 2009–2019, and all the changes (i.e., death and new in-growth) during the periods were recorded. We found the time stabilities of survival functions between periods in advance. The results showed that the parametric survival analysis with Weibull distribution successfully yielded the mortality rate, mortality probability, and survival probability in each plot. Finally, we predicted the future age class distribution of living and dead trees of each plot based on the survival analyses results and discussed its management implications. We recommended that the estimated mean lifetime was facilitated to make decisions on the selection of harvesting trees in the uneven forest management based on selective cutting.


Introduction
Changes in the survival probability of forests will have environmental [1] and severe economic consequences.Survival probability changes can influence forest management decisions [2,3].There is a strong need for long-term tree mortality studies [4][5][6][7], considering the threat of global climate change, which is forecasted to increase tree mortality in forests worldwide [8].
Because of the difficulty in detecting the exact ages of trees, researchers have mainly used the Diameter at Breast Height (DBH), dominant height, basal area, growth rate, or competition index, avoiding age-based measures, to determine mortality probabilities.Generally, a tree's full survival time is unknown; unlike in logistic regressions, the right censoring and left truncation approaches of survival analysis can be applied to handle these data.For the first time, [9] suggested survival analysis techniques to solve forest mortalities, even though these applications were limited to even-aged forest stands  [10,11].Woodall [12] carried out a survival analysis based on the Diameter at Breast Height (DBH).Finally, he emphasized the importance of the exact age on the survival analysis of forests.
The principles of the approach presented in this paper are related to the basic idea shown by Hiroshima [13], who applied more flexible techniques of survival analysis based on non-parametric and parametric approaches.This method is based on annual tree ring measurements with a semi-nondestructive device as well as periodic inventory data of secondary natural forests in Japan.However, except for this, studies focusing on the agebased survival analysis of unevenly aged natural forests are comparatively rare because they require special methods to detect trees' ages and make long-term and large-scale observations as well.In this study, we used a particular device to detect tree age and take advantage of the long-term research site in the preserved area, located in the University of Tokyo Hokkaido Forest (UTHF), profoundly affected by wind-induced forest mortality.Overall, the purpose of this study was to perform a survival analysis of individual tree populations in the natural forest stands of the UTHF, Japan, to predict the age distributions of living and dead trees.First, we applied survival analysis techniques to each plot to make non-parametric estimations.Secondly, parametric Weibull estimates were developed to predict future age class distributions.Specifically, we showed how the stand age estimations differ with mean lifetime after ingrowth in each plot.

Study Area
The study was carried in the UTHF, which is located in Furano (central Hokkaido, northern Japan).
In this study, we chose two permanent plots (0.25 ha each), which are located at elevations of 681 and 580 m in the preserved area (Figure 1).These middle elevations represent the range of the species composition as well as stand structure, with a wind-damaged history.The typical vegetation types of these plots are a coniferous and broad-leaved mixed forest dominated by Abies sachalinensis and Acer species.

Collection of Tree Age Data
Across the two plots, we detected the numbers of tree rings at breast height (1.3 m) of 451 trees in 2019 using a semi-nondestructive device called the RESISTOGRAPH, developed by [14].The 10-year observation period was set as 2009-2019.All the target trees were alive and had DBH ≥ 5 cm in 2009, and some of these trees had died by 2019.The live/dead statuses and new ingrowth trees were counted at the end of each observation period.The field measurement data of the tree rings in 2019 were extracted via the DE-COM software to detect annual rings.Figure 2 describes the annual rings detected in one of the alive Abies trees using the DECOM software.The RESOSTOGRAPH measurements were ineffective for badly rotten and center-decayed trees.Cross-sectional wood discs were cut in the field to count the annual rings.We manually counted the annual rings for those species and estimated the equation of the regression between the radius of the wood disk (x) and the number of annual rings (y) [15].The age after ingrowth was determined by inputting the radius data into simple three-dimensional equations fitted to a scatter diagram.This method was followed for 25% of the sample trees.
We analyzed/estimated 237 and 214 trees of plot numbers 5225 and 5240, respectively.The relevant age at the beginning of the observation periods in 2009 was determined by deducting ten years from the deduced tree age in 2019 of each tree.Finally, the tree ages were classified with ten years per age class.

Survival Analysis of Target Trees
The survival probability functions were the essential elements of this analysis, as they reflect the overall performance of the two plots.These functions were developed based on previous studies [16][17][18].The survival data demonstrate that there were both dead and alive trees in the study plots at the end of the observation period.In these instances, the survival time is described as both "censored" and "truncated".Right censoring happens when a target leaves the study before an event occurs, or the follow-up ends before it happens.On the other hand, left truncation happens when an object is not observed from the start of the study but instead enters the observation period in this study later.
The Kaplan-Meier survival estimator is an important tool for analyzing censored data, and the Kaplan-Meier curve [19] can be used to estimate the distribution of tree mortality.
The survival analysis was performed to describe the survival probability distribution by following the above studies.The survival function of the Kaplan-Meier (KM) estimate is as follows: where t r ∧ represents the survival function and represents the hazard function for the t th age class.
Besides non-parametric analysis, the Weibull distribution was applied to predict the future age distributions of living and dead trees.This parametric analysis estimates the parameters m and k of the Weibull distribution to check the fitness of Kaplan-Meier estimates.(2) Equations ( 3) and (4) were used to predict the living and dead tree distributions after the second age class of each plot.
where bt and ct are the numbers of living trees and dead trees in the current period by age class t based on the number of living trees at, dead trees dt and corresponding tree mortality rate pt derived from Weibull analysis for the previous period.

Results
The species compositions of the living and dead trees of the two plots were different, so we avoided combining the tree data in the two plots to carry out the following survival analyses.We found that these Kaplan-Meier curves were time stable in both plots through the log-rank test [15]; therefore, these survival probabilities were used to carry out the future predictions in the following parametric analyses.
The Weibull distribution was applied to smooth the stepwise values of the non-parametric analysis for prediction purposes.Table 1 shows the estimated parameters of m and k and the means and standard deviations of the Weibull distributions of two plots.Higher mean and standard deviation values were found for plot 5225.The mean value represents the mean lifetime after ingrowth, and according to the Weibull estimates, the mean lifetimes of 5225 and 5240 are 137 and 95 years, respectively.For validation purposes, the living and dead tree distributions of each plot in the observation period were calculated based on Equations ( 3) and ( 4).The studied plots were characterized by uneven age class distributions of trees.However, the shape of the age class distribution varied between the two plots.The proportion of living and dead trees of 5225 declined gradually toward the older age classes, while 5240 represented a slight difference in the age class frequencies for the younger age classes.For t ≥ 2, the total numbers of Weibull calculations were 157 (plot 5225) and 159 (plot 5240), while the observations were 155 (plot 5225) and 158 (plot 5240), respectively.The estimation of surviving trees ended up showing 1.29 and 0.63% error ratios for 5225 and 5240 and represented the characteristics of the age distributions correctly.The error ratios for the dead trees in 5225 and 5240 were 16.67 and 0%, respectively.The error ratios based on the parametric analysis nearly agreed with the observed data, which shows that the validation results were sufficiently precise.
Based on the estimated pt during 2009-2019 and assuming that these mortality rates were stable over time, the age class distribution in the next period, 2019-2029, was predicted for each plot.The resulting living and dead tree distributions are shown in Figure 3.The predicted numbers of living trees were less than in the 2009-2019 period, especially in the second, third, and fourth age classes.By contrast, the living tree projections for the second and third age classes were higher than those observed in 5240.

Discussion
The estimates of the shape parameters differed in the range of k < 1 for 5225 and k > 1 for 5240, which led to the differences between the shapes of the curves in the two plots.Because the dead trees in 5240 tended to concentrate in the older age classes compared with 5225, the estimates of k for 5240 showed higher values than 5225.These differences in shapes resulted in a difference in mean lifetime of 137 years in 5225 and 95 years in 5240, which was mainly caused by the difference in mortality rates between the plots.The older tree death of 5240 led to a lower mean lifetime than in 5225, which also shows that a higher number of trees could survive in 5225 than 5240.The weighted means of these biological lifetimes were calculated as 171 and 173 years in two plots based on the species and number of trees.Thus, the above mean lifetimes of 137 and 95 years and also the maximum stand ages of 159 and 149 years in both plots were relatively near to the mean biological lifetimes of 171 and 173 years.These facts imply the matured states of these stands, which led to the stable states of the survival probabilities over time and enhance the future predictions based on the estimated survival probabilities.
Forest managers can rely on these predicted age class distributions in decision-making, such as harvesting tree selection based on the Stand-based Silvicultural Management System (SSMS).The resulting future prediction of dead trees concerning age classes will facilitate, as an indicator to find the suitable trees for the harvesting.The productivity of harvested trees can be enhanced by this type of early identification of likely-to-die or decaying trees.
Thus, future predictions of the age class distribution of each plot such as that in this study can be used in site-specific management plans to identify living and dead trees in each age class, though it is essential to carry out dendrological survival analysis instead of whole-stand analysis.

Conclusions
In this study, we succeeded in estimating survival probabilities in both nonparametric and parametric ways.We also showed the feasibility of future prediction using parametric estimates of survival probabilities.The results imply that natural forest management based on selective cutting such as SSMS in UTHF could be practiced more efficiently by considering the estimated mean lifetimes derived from survival analysis.
However, the survival probabilities estimated in this study should be used carefully for long-term predictions of forest dynamics because they do not include the effects of catastrophic disturbances, which can often have significant influences on forests.Therefore, in our future work, it will be essential to incorporate these variables to enhance the practical applicability of survival probability models.

Citation:
Wijenayake, P.R.; Hiroshima, T. Prediction of Tree Age Distribution Based on Survival Analysis in Natural Forests: A Case Study of Preserved Permanent Plots in the University of Tokyo Hokkaido Forest, Northern Japan.Environ.Sci.Proc.2021, 3, 50.https://doi.org/10.3390/IECF2020-08077 Academic Editors: Angela Lo Monaco, Cate Macinnis-Ng and Om P. Rajora Published: 13 November 2020 Publisher's Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.Copyright: © 2020 by the authors.Licensee MDPI, Basel, Switzerland.This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Figure 1 .
Figure 1.Maps showing the locations of study area and plots. (a) University of Tokyo Hokkaido Forest (UTHF) of the Hokkaido island of Japan, (b) preserved area and permanent plots of the UTHF, and (c) contour map of the two preserved permanent plots of UTHF.

Figure 2 .
Figure 2. Manually detected annual rings of Abies sachalinensis tree of plot 5240, which were counted after detection with the DECOM software.

Table 1 .
Weibull probability distributions for each plot.
m = scale parameter, k = shape parameter.