Next Article in Journal
Sustainable Forest Management Evaluation Using Carbon Credits: From Production to Environmental Forests
Next Article in Special Issue
Bark Thickness and Heights of the Bark Transition Area of Scots Pine
Previous Article in Journal
Assessment of Tree Diameter Estimation Methods from Mobile Laser Scanning in a Historic Garden
Previous Article in Special Issue
Soil Temperature in Disturbed Ecosystems of Central Siberia: Remote Sensing Data and Numerical Simulation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Department of Forest Science, Graduate School of Agricultural and Life Sciences, The University of Tokyo, Bunkyo-ku, Tokyo 113-8657, Japan
2
Department of Global Agricultural Sciences, Graduate School of Agricultural and Life Sciences, The University of Tokyo, Bunkyo-ku, Tokyo 113-8657, Japan
*
Author to whom correspondence should be addressed.
Submission received: 30 June 2021 / Revised: 24 July 2021 / Accepted: 26 July 2021 / Published: 30 July 2021

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). 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 DECOMTM 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 DECOMTM 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). This method was followed for 20% of the sample trees. Table 3 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):
Pr ( t 1 < T t T > t 1 ) = p t
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:
Pr ( T > t T > t 1 ) = 1 p t .
Following previous methods [45,46,47] and considering Equations (1) and (2), we described the likelihood function (L) of the observation as follows:
L = t Pr ( t 1 < T t T > t 1 ) d t Pr ( T > t T > t 1 ) a t = t p t d t 1 p t a t ,
where dt is the number of dead trees and at 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:
Pr ( t 1 < T t ) = q t .
The survival probability ( r t ) in the t-th age class was defined as follows:
Pr ( T > t 1 ) = r t .
Therefore, based on Equation (1), the mortality rate ( p t ) can be expressed as follows:
p t = q t r t .
The maximum likelihood estimators of p t can be calculated by the first-order derivation of Equation (3), as shown in Equation (7):
p t = d t a t + d t .
Considering that Equations (5) and (6) can be combined as follows:
r t r t + 1 = q t ,
the survival function can be converted into:
r t + 1 = r t q t = r t 1 q t r t = r t 1 p t .
Equation (9) can also be expressed as follows:
r t = k < t 1 p k .
The survival function is described by the following formula:
r t = k < t 1 p k = k < t 1 d k a k + d k ,
where r t represents the survival function and   d k / ( a k + d k ) represents the hazard function in the t-th age class.
This 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:
Log rank statistic = O E 2 / Var   O E .
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:
Wilcoxon   test   statistics = f w ( t ( f ) ) ( O E ) 2 j w ( t ( f ) ) ( O E ) .
It weights the observed minus expected score at the time tf by the number at risk, w(tf) overall groups at time tf.
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:
Mean   lifetime = t q t / q t .
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. 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), 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 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 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). 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). 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.

Author Contributions

Data analysis, writing—original draft preparation, P.R.W.; Supervision, review and editing, T.H. Both authors have read and agreed to the published version of the manuscript.

Funding

This study was funded by the JSPS KAKENHI, Grant Number JP19K06142.

Data Availability Statement

The datasets relevant to the current study are available from the corresponding author on reasonable request.

Acknowledgments

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

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Seidl, R.; Schelhaas, M.J.; Rammer, W.; Verkerk, P.J. Increasing forest disturbances in Europe and their impact on carbon storage. Nat. Clim. Chang. 2014, 4, 806–810. [Google Scholar] [CrossRef] [Green Version]
  2. Neuner, S.; Knoke, T. Economic consequences of altered survival of mixed or pure Norway spruce under a dryer and warmer climate. Clim. Chang. 2017, 140, 519–531. [Google Scholar] [CrossRef]
  3. Griess, V.C.; Knoke, T. Bioeconomic modeling of mixed Norway spruce-European beech stands: Economic consequences of considering ecological effects. Eur. J. For. Res. 2013, 132, 511–552. [Google Scholar] [CrossRef]
  4. Paul, C.; Brandl, S.; Friedrich, S.; Falk, W.; Härtl, F.; Knoke, T. Climate change and mixed forests: How do altered survival probabilities impact economically desirable species proportions of Norway spruce and European beech? Ann. For. Sci. 2019, 76, 14. [Google Scholar] [CrossRef] [Green Version]
  5. Adams, H.D.; MacAlady, A.K.; Breshears, D.D.; Allen, C.D.; Stephenson, N.L.; Saleska, S.R.; Huxman, T.E.; McDowell, N.G. Climate-induced tree mortality: Earth system consequences. EOS Trans. Am. Geophys. Union 2010, 91, 153–154. [Google Scholar] [CrossRef]
  6. Dietze, M.C.; Moorcroft, P.R. Tree mortality in the eastern and central United States: Patterns and drivers. Glob. Chang. Biol. 2011, 17, 3312–3326. [Google Scholar] [CrossRef]
  7. Pfeifer, E.M.; Hicke, J.A.; Meddens, A.J.H. Observations and modeling of aboveground tree carbon stocks and fluxes following a bark beetle outbreak in the western United States. Glob. Chang. Biol. 2011, 17, 339–350. [Google Scholar] [CrossRef]
  8. Trumbore, S.; Brando, P.; Hartmann, H. Forest health and global change. Science 2015, 349, 814–818. [Google Scholar] [CrossRef] [Green Version]
  9. Runkle, J.R. Canopy tree turnover in old-growth mesic forests of eastern north America. Ecology 2000, 81, 554–567. [Google Scholar] [CrossRef]
  10. Harcombe, P.A.; Bill, C.J.; Fulton, M.; Glitzenstein, J.S.; Marks, P.L.; Elsik, I.S. Stand dynamics over 18 years in a southern mixed hardwood forest, Texas, USA. J. Ecol. 2002, 90, 947–957. [Google Scholar] [CrossRef]
  11. Van Mantgem, P.J.; Stephenson, N.L.; Byrne, J.C.; Daniels, L.D.; Franklin, J.F.; Fulé, P.Z.; Harmon, M.E.; Larson, A.J.; Smith, J.M.; Taylor, A.H.; et al. Widespread increase of tree mortality rates in the Western United States. Science 2009, 323, 521–524. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Allen, C.D.; Macalady, A.K.; Chenchouni, H.; Bachelet, D.; McDowell, N.; Vennetier, M.; Kitzberger, T.; Rigling, A.; Breshears, D.D.; Hogg, E.H.; et al. A global overview of drought and heat-induced tree mortality reveals emerging climate change risks for forests. For. Ecol. Manag. 2010, 259, 660–684. [Google Scholar] [CrossRef] [Green Version]
  13. McDowell, N.G.; Allen, C.D. Darcy’s law predicts widespread forest mortality under climate warming. Nat. Clim. Chang. 2015, 5, 669–672. [Google Scholar] [CrossRef]
  14. Woodall, C.W.; Grambsch, P.L.; Thomas, W. Applying survival analysis to a large-scale forest inventory for assessment of tree mortality in Minnesota. Ecol. Model. 2005, 189, 199–208. [Google Scholar] [CrossRef]
  15. Monserud, R. Simulation of Forest Tree Mortality. For. Sci. 1976, 22, 438–444. [Google Scholar] [CrossRef]
  16. Hamilton, D.A. A logistic model of mortality in thinned and unthinned mixed conifer stands of northern Idaho. For. Sci. 1986, 32, 989–1000. [Google Scholar] [CrossRef]
  17. Kobe, R.K.; Coates, K.D. Models of sapling mortality as a function of growth to characterize interspecific variation in shade tolerance of eight tree species of northwestern British Columbia. Can. J. For. Res. 1997, 27, 227–236. [Google Scholar] [CrossRef]
  18. Buford, M.A.; Hafley, W.L. Probability distributions as models for mortality. For. Sci. 1985, 31, 331–341. [Google Scholar]
  19. Finney, D.J. Probit Analysis, 3rd ed.; Cambridge University Press: London, UK, 1971. [Google Scholar]
  20. Kalhfleisch, J.G.; Prentice, R.L. The Statistical Analysis of Failure Time Data; John Wiley: New York, NY, USA, 1980. [Google Scholar]
  21. Fujikake, I. An Analyses of Harvesting Activities in a Forest Management Entity Using Harvesting Age Distribution (Tentative Translation by Author). Ph.D. Thesis, Kyoto University, Kyoto, Japan, 2000. [Google Scholar]
  22. Kleinbaum, D.G.D.; Klein, M. Survival Analysis: A Self-Learning Text, 3rd ed.; Statistics for Biology and Health; Springer: Berlin/Heidelberg, Germany, 2011. [Google Scholar]
  23. Kaplan, E.L.; Meier, P. Nonparametric Estimation from Incomplete Observations. J. Am. Stat. Assoc. 1958, 53, 457–481. [Google Scholar] [CrossRef]
  24. Waters, W.E. Life-table approach to analysis of insect impact. J. For. 1969, 67, 300–304. [Google Scholar]
  25. Morse, B.W.; Kulman, H.M. Plantation white spruce mortality: Estimates based on aerial photography and analysis using a life-table format. Can. J. For. Res. 1984, 14, 195–200. [Google Scholar] [CrossRef]
  26. Zhang, S.; Amateis, R.L.; Burkhart, H.E. Constraining individual tree diameter increment and survival models for loblolly pine plantations. For. Sci. 1997, 43, 413–423. [Google Scholar] [CrossRef]
  27. Wyckoff, P.H.; Clark, J.S. Predicting tree mortality from diameter growth: A comparison of maximum likelihood and Bayesian approaches. Can. J. For. Res. 2000, 30, 156–167. [Google Scholar] [CrossRef]
  28. Nothdurft, A. Spatio-temporal prediction of tree mortality based on long-term sample plots, climate change scenarios and parametric frailty modeling. For. Ecol. Manag. 2013, 291, 43–54. [Google Scholar] [CrossRef]
  29. Neuner, S.; Albrecht, A.; Cullmann, D.; Engels, F.; Griess, V.C.; Hahn, W.A.; Hanewinkel, M.; Härtl, F.; Kölling, C.; Staupendahl, K.; et al. Survival of Norway spruce remains higher in mixed stands under a dryer and warmer climate. Glob. Chang. Biol. 2015, 21, 935–946. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  30. Neumann, M.; Mues, V.; Moreno, A.; Hasenauer, H.; Seidl, R. Climate variability drives recent tree mortality in Europe. Glob. Chang. Biol. 2017, 23, 4788–4797. [Google Scholar] [CrossRef]
  31. Brandl, S.; Paul, C.; Knoke, T.; Falk, W. The influence of climate and management on survival probability for Germany’s most important tree species. For. Ecol. Manag. 2020, 458, 117652. [Google Scholar] [CrossRef]
  32. Hiroshima, T. Applying age-based mortality analysis to a natural forest stand in Japan. J. For. Res. 2014, 19, 379–387. [Google Scholar] [CrossRef]
  33. Manion, P.D.; Griffin, D.H. Large landscape scale analysis of tree death in the Adirondack Park, New York. For. Sci. 2001, 47, 542–549. [Google Scholar] [CrossRef]
  34. Zens, M.S.; Peart, D.R. Dealing with death data: Individual hazards, mortality and bias. Trends Ecol. Evol. 2003, 18, 366–373. [Google Scholar] [CrossRef]
  35. The University of Tokyo Hokkaido Forest. 2019. Available online: http://www.uf.a.u-tokyo.ac.jp/files/gaiyo_hokkaido.pdf (accessed on 26 July 2020).
  36. Asahi, M. Studies on the classification of forest soils in the Tokyo University Forest, Hokkaido. Bull. Tokyo Univ. For. 1963, 58, 1–132. [Google Scholar]
  37. Watanabe, S.; Sasaki, S. The silvicultural management system in temperate and boreal forests—A case history of the Hokkaido Tokyo University Forest. Can. J. For. Res. 1994, 24, 1176–1185. [Google Scholar] [CrossRef]
  38. Owari, T.; Matsui, M.; Inukai, H.; Kaji, M. Stand structure and geographic conditions of natural selection orests in central Hokkaido, Northern Japan. J. For. Plan. 2011, 16, 207–214. [Google Scholar]
  39. Rinn, F. Eine neue Bohrmethode zur Holzuntersuchung. Holz-Zentralblatt 1989, 34, 529–530. [Google Scholar]
  40. Wijenayake, P.R.; Hiroshima, T. Survival analyses of individual tree populations in natural forest stands to evaluate the maturity of forest stands: A case study of preserved forests in Northern Japan ). J. For. Plan. 2021. [Google Scholar] [CrossRef]
  41. Cox, D.R.; Oakes, D. Analysis of Survival Data; Chapman and Hall: London, UK, 1984. [Google Scholar]
  42. Crowder, M.J.; Kimber, A.C.; Smith, R.L. Statistical Analysis of Reliability Data.; Chapman and Hall: London, UK, 1994; Volume 27. [Google Scholar]
  43. Klein, J.P.; Moeschberger, M.L. Survival Analysis: Techniques for Censored and Truncated Data; Springer: New York, NY, USA, 1997. [Google Scholar]
  44. Cox, D.R. Regression Models and Life-Tables. J. R. Stat. Soc. 1972, 34, 187–220. [Google Scholar] [CrossRef]
  45. Fujikake, I. Estimation of Gentan probability based on the forest resource table. Proc. Inst. J. Stat. Math. 2003, 51, 95–109. [Google Scholar] [CrossRef]
  46. Hiroshima, T. Study on the methodology of calculating mean and variance of felling age in forest planning. Jpn. J. For. Plan. 2006, 40, 139–149. [Google Scholar] [CrossRef]
  47. Tiryana, T.; Tatsuhara, S.; Shiraishi, N. Empirical models for estimating the stand biomass of teak plantations in Java, Indonesia (Multipurpose Forest Management). J. For. Plan. 2011, 16, 35–44. [Google Scholar] [CrossRef]
  48. Mantel, N. Evaluation of survival data and two new rank order statistics arising in its consideration. Cancer Chemother. Rep. 1966, 50, 163–170. [Google Scholar]
  49. Gehan, E. A generalised Wilcoxon test for comparing arbitrarily singly-censored samples. Biometrika 1965, 52, 203–223. [Google Scholar] [CrossRef]
  50. Hiura, T.; Go, S.; Iijima, H. Long-term forest dynamics in response to climate change in northern mixed forests in Japan: A 38-year individual-based approach. For. Ecol. Manag. 2019, 449, 117469. [Google Scholar] [CrossRef]
  51. Sato, T. Time trends for genetic parameters in progeny test of Abies sachalinensis (Fr. Schm.) Mast. Silvae Genet. 1994, 43, 304–307. [Google Scholar]
  52. Hiura, T.; Sano, J.; Konno, Y. Age structure and response to fine-scale disturbances of Abies sachalinensis, Picea jezoensis, Picea glehnii, and Betula ermanii growing under the influence of a dwarf bamboo understory in northern Japan. Can. J. For. Res. 1996, 26, 289–297. [Google Scholar] [CrossRef]
  53. Nakashizuka, T. Population dynamics of coniferous and broad-leaved trees in a Japanese temperate mixed forest. J. Veg. Sci. 1991, 2, 413–418. [Google Scholar] [CrossRef]
  54. Selvin, S. Statistical Analysis of Epidemiologic Data; Oxford University Press: Oxford, UK, 2004. [Google Scholar]
  55. Yamamoto, H. Natural forest management based on selection cutting and natural regeneration. In Proceedings of the IUFRO International Workshop on sustainable Forest Managements, Furano, Japan, 17–21 October 1994; pp. 10–22. [Google Scholar]
  56. Takahashi, N. The Stand-Based Silvicultural Management System: Its Theory and Practices, 1st ed.; Log Bee: Sapporo, Japan, 2001. [Google Scholar]
  57. Franklin, J.F.; Spies, T.A.; Van Pelt, R.; Carey, A.B.; Thornburgh, D.A.; Berg, D.R.; Lindenmayer, D.B.; Harmon, M.E.; Keeton, W.S.; Shaw, D.C.; et al. Disturbances and structural development of natural forest ecosystems with silvicultural implications, using Douglas-fir forests as an example. For. Ecol. Manag. 2002, 155, 399–423. [Google Scholar] [CrossRef]
  58. Daniel, T.W.; Helms, J.A.; Baker, S.F. Principles of Silviculture, 2nd ed.; McGraw-Hill: New York, NY, USA, 1979. [Google Scholar]
Figure 1. Locations of the study area and plots. (a) University of Tokyo Hokkaido Forest (UTHF) of Hokkaido island of Japan; (b) preserved area and permanent plots of the UTHF; (c) contour map of the four preserved permanent plots of UTHF.
Figure 1. Locations of the study area and plots. (a) University of Tokyo Hokkaido Forest (UTHF) of Hokkaido island of Japan; (b) preserved area and permanent plots of the UTHF; (c) contour map of the four preserved permanent plots of UTHF.
Forests 12 01014 g001
Figure 2. Manually detected annual rings of A. sachalinensis tree in plot 5240 detecting using DECOM TM software.
Figure 2. Manually detected annual rings of A. sachalinensis tree in plot 5240 detecting using DECOM TM software.
Forests 12 01014 g002
Figure 3. Age class distributions of living and dead trees of major species: (a) A. sachalinensis, (b) P. jezoensis, and (c) all broad-leaved species.
Figure 3. Age class distributions of living and dead trees of major species: (a) A. sachalinensis, (b) P. jezoensis, and (c) all broad-leaved species.
Forests 12 01014 g003
Figure 4. Kaplan–Meier curves showing tree survival probabilities during the three observational periods: (a) A. sachalinensis, (b) P. jezoensis, and (c) all broad-leaved species.
Figure 4. Kaplan–Meier curves showing tree survival probabilities during the three observational periods: (a) A. sachalinensis, (b) P. jezoensis, and (c) all broad-leaved species.
Forests 12 01014 g004
Figure 5. Kaplan–Meier curves showing tree survival probabilities of A. sachalinensis and P. jezoensis during the three observational periods: (a) period 1, (b) period 2, and (c) period 3.
Figure 5. Kaplan–Meier curves showing tree survival probabilities of A. sachalinensis and P. jezoensis during the three observational periods: (a) period 1, (b) period 2, and (c) period 3.
Forests 12 01014 g005
Figure 6. Kaplan–Meier curves for each period for conifers and broad-leaved trees: (a) period 1, (b) period 2, and (c) period 3 (i and ii denotes the starting age classes of “practical” mean lifetime of each broad-leaved and coniferous species, respectively).
Figure 6. Kaplan–Meier curves for each period for conifers and broad-leaved trees: (a) period 1, (b) period 2, and (c) period 3 (i and ii denotes the starting age classes of “practical” mean lifetime of each broad-leaved and coniferous species, respectively).
Forests 12 01014 g006
Table 1. Stand characteristics of selected plots.
Table 1. Stand characteristics of selected plots.
PlotPlot Size (ha)Elevation (m)No. of Target Trees in Period 3Slope AspectMean Slope Angle (°)
52030.40580327Southwest15
52240.25570150Southwest18
52250.25690237Southwest20
52400.25600214Southwest25
Table 2. Regression models developed based on RESISTOGRAPH measurements.
Table 2. Regression models developed based on RESISTOGRAPH measurements.
Tree SpeciesNo. of Tree SamplesEquationCoefficient of Determination
A. sachalinensis51y = 0.0000008w3 + 0.0005w2 + 0.1995w + 32.2940.8379
P. jezoensis59y = 0.000002w3 + 0.0008w2 + 0.2545w + 30.320.8506
A. ukurunduense41y = 0.00001w3 + 0.0043w2 + 0.1416w + 41.9910.7668
T. japonica32y = 0.000002w3 − 0.0012w2 + 0.5815w + 13.1890.7803
Other broad-leaved species45y = 0.000005w3 + 0.0023 w2 + 0.0088 w + 42.3860.8592
Table 3. Methodologies used to determine the age of each tree species.
Table 3. Methodologies used to determine the age of each tree species.
SpeciesCensus DataDirect
AResistograph Measurements
Developed
ARegression Models
Total No. of Trees
Living TreesDead TreesLiving TreesDead TreesLiving TreesDead TreesLiving TreesDead Trees
A. sachalinensis113958164816 219 41
P. jezoensis1845511496 117 21
Broad-leaved trees2896776215510 423 98
Percentage45.708.7120.575.2216.543.48
Table 4. Summary of the dominant tree species in the four plots in the three observation periods.
Table 4. Summary of the dominant tree species in the four plots in the three observation periods.
SpeciesNo. of Trees (%)
Period 1Period 2Period 3
Living TreesDead TreesLiving TreesDead TreesLiving TreesDead Trees
Conifer
A. sachalinensis201 (21.90)22 (2.40)215 (22.16)22 (2.27)219 (23.60)41 (4.42)
P. jezoensis148 (16.12)14 (1.53)135 (13.92)16 (1.65)117 (12.61)22 (2.37)
Other conifers4 (0.44)2 (0.22)1 (0.10)3 (0.31)1 (0.11)0 (0.00)
Total conifer353 (38.45)38 (4.14)351 (36.19)41 (4.23)337 (36.31)63 (6.79)
Broad-leaved
Ulmus laciniata56 (6.10)8 (0.87)52 (5.36)10 (1.03)46 (4.96)8 (0.86)
Sorbus commixta35 (3.81)3 (0.33)27 (2.78)11 (1.13)25 (2.69)8 (0.86)
Acer ukurunduense21 (2.29)9 (0.98)48 (4.95)11 (1.13)50 (5.39)19 (2.05)
Acer mono var. myrii56 (6.10)2 (0.22)57 (5.88)6 (0.62)54 (5.82)7 (0.75)
Tilia japonica174 (18.95)29 (3.16)171 (17.63)52 (5.36)153 (16.49)32 (3.45)
Other broad-leaved115 (12.53)19 (2.07)113 (11.65)20 (2.06)101 (10.88)25 (2.69)
Total broad-leaved457 (49.78)70 (7.63)468 (48.25)110 (11.34)429 (46.23)99 (10.67)
Total810 (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.
Table 5. Log-rank test results comparing the three observational periods.
Table 5. Log-rank test results comparing the three observational periods.
PeriodTree Species
A. sachalinensisP. jezoensisAll Brzoad-Leaved Species
Log-Rank TestWilcoxon TestLog-Rank TestWilcoxon TestLog-Rank TestWilcoxon Test
p-ValueChi-Square Valuep-ValueChi-Square Valuep-ValueChi-Square Valuep-ValueChi-Square Valuep-ValueChi-Square Valuep-ValueChi-Square Value
Periods 1 and 20.91610.01110.65420.20070.50230.45020.36740.81240.0194 *5.46370.07723.1236
Periods 2 and 30.10342.65240.56850.32520.15382.03370.25491.29610.26911.22150.06103.5101
* Indicates a significant difference (p < 0.05).
Table 6. Log-rank test results of Abies sachalinensis and Picea jezoensis in each period.
Table 6. Log-rank test results of Abies sachalinensis and Picea jezoensis in each period.
A. sachalinensis and P. jezoensis
Period 1Period 2Period 3
Statistical Testp-ValueChi-Square Valuep-ValueChi-Square Valuep-ValueChi-Square Value
Log-rank test0.49850.45810.57510.31420.99590.0000
Wilcoxon test0.48270.49280.57200.31930.86920.0271
Table 7. Log-rank test and Wilcoxon test comparing the conifer and broad-leaved trees in each period.
Table 7. Log-rank test and Wilcoxon test comparing the conifer and broad-leaved trees in each period.
Conifer and broad-Leaved Species
Period 1Period 2Period 3
Statistical Tzstp-ValueChi-Square Valuep-ValueChi-Square Valuep-ValueChi-Square Value
Log-rank test0.16141.96070.0097 *6.68700.0381 *4.3018
Wilcoxon test0.0386 *4.27970.0463 *3.97190.30351.0586
* indicates a significant difference (p < 0.05).
Table 8. Mean lifetime and standard deviation (in parenthesis) of two groups of trees (coniferous and broad-leaves) in each period.
Table 8. Mean lifetime and standard deviation (in parenthesis) of two groups of trees (coniferous and broad-leaves) in each period.
PeriodConifer SpeciesBroad-Leaved Species
Period 139 (4)44 (5)
Period 236 (4)36 (5)
Period 343 (5)42 (5)
Table 9. The predicted practical mean lifetime and corresponding diameter at breast height value of period 3.
Table 9. The predicted practical mean lifetime and corresponding diameter at breast height value of period 3.
SpeciesPractical Mean Lifetime in Years (SD)Practical Mean DBH in cm (SD)
Conifer species72 (7)31 (7)
Broad-leaved species80 (5)36 (8)
Diameter at breast = DBH; standard deviation = SD.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Wijenayake, P.R.; Hiroshima, T. Age-Based Survival Analysis of Coniferous and Broad-Leaved Trees: A Case Study of Preserved Forests in Northern Japan. Forests 2021, 12, 1014. https://0-doi-org.brum.beds.ac.uk/10.3390/f12081014

AMA Style

Wijenayake PR, Hiroshima T. Age-Based Survival Analysis of Coniferous and Broad-Leaved Trees: A Case Study of Preserved Forests in Northern Japan. Forests. 2021; 12(8):1014. https://0-doi-org.brum.beds.ac.uk/10.3390/f12081014

Chicago/Turabian Style

Wijenayake, Pavithra Rangani, and Takuya Hiroshima. 2021. "Age-Based Survival Analysis of Coniferous and Broad-Leaved Trees: A Case Study of Preserved Forests in Northern Japan" Forests 12, no. 8: 1014. https://0-doi-org.brum.beds.ac.uk/10.3390/f12081014

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop