Climatic information improves statistical individual-tree mortality models for three key species of Sichuan Province, China

Abstract

• Key message

Climate variables improve individual-tree mortality models for fir, oak and birch.

• Context

Climate is considered as an important driver of tree mortality, but few studies have included climate factors in models to explore their importance for modelling individual-tree mortality.

• Aims

To measure the performance of climate-based models, we built individual-tree mortality models using individual, stand, and climate variables for fir (Abies faxoniana Rehd. et Wils.), oak (Quercus aquifolioides Rehd. et Wils.), and birch (Betula albo-sinensis Burk.) in Southwest China, and explored the corresponding effects on tree death.

• Methods

We developed tree mortality models based on 287 permanent plots from the Sichuan Forest Inventory data, and compared the models based on variables of individual (I), stand (S), and climate (C) levels, and different combinations (I + S, I + C, S + C, I + S + C) among these groups to improve model performance. We employed relative Akaike information criterion (AIC), area under receiver operating characteristic curve (AUC), and Hosmer-Lemeshow’s goodness-of-fit statistic for model evaluation and validation.

• Results

We found that tree mortalities of the three species could be better predicted (AUC > 0.8) by carefully selecting variables at three ecological scales (individual, stand, and regional climate). Our results suggest that the higher mortality of the object trees occurs when they endure lower radial growth of the previous years, more intensive competition, and moderate canopy cover (for birch), while lower mortality was seen in an appropriate range of climate conditions and at higher stand canopy cover (for fir and oak).

• Conclusion

The results have significance for incorporating the effects of a changing climate into mortality models.

1 Introduction

Tree death is a process that removes individuals from a stand, providing additional resources back to the ecosystem through the decomposition of dead trees. Thus, tree mortality plays a key role in shaping forest structure and stands dynamics, and drives a number of ecosystem functions and processes, such as carbon sequestrations and emissions, nutrient cycling, and water production (Franklin et al. 1987). Tree mortality can often result from catastrophic events, such as wildfire, hurricane, snow storm, insect outbreaks, and geologic hazards (Van Mantgem et al. 2009; Peng et al. 2011; Brando et al. 2014), or non-catastrophic events related to resource competitions and growth conditions, such as light regime, climate, water availability, soil nutrients, and senescence (Fridman and Ståhl 2001; Crecente-Campo et al. 2009). Although non-catastrophic tree death may have little effect on stand productivity, as well demonstrated by the thinning-response hypothesis especially for even-aged forests (Skovsgaard and Vanclay 2008), dead trees have a large impact on litter production and soil carbon storage, thus significantly affecting ecosystem respiration and net ecosystem production (Clark and Clark 1996).

Tree mortality is also indispensable for assessing forest health status and managing forest resources sustainably (Wunder et al. 2007). Modeling tree mortality has been a challenge due to the complicated interactions among biotic and abiotic factors relating to tree growth and death, and our insufficient knowledge of the processes of tree aging and death. Non-catastrophic tree mortality (also called natural mortality) resulting from competition for light, water and soil nutrients within a stand is widely observed in various forests, while catastrophic tree mortality is more abrupt and often occurs at more local scales (Peet and Christensen 1987; Crecente-Campo et al. 2009). Therefore, most tree-mortality models have focused on non-catastrophic mortality because of the highly variable and stochastic nature of the catastrophic events (Fridman and Ståhl 2001; Bigler and Bugmann 2004; Crecente-Campo et al. 2009). Individual-tree-based mortality models not only show details of each tree for improving forest management but also provide information about the processes and factors related to individual tree death (Monserud and Sterba 1999; Wunder et al. 2008; Adame et al. 2010). Moreover, many ecosystem models have not explicitly considered tree mortality due to the lack of information on tree death and its relationships with climate factors (Wang et al. 2012).

Biotic and abiotic factors may contribute to tree mortality at different ecological scales, such as individual susceptibility, stand features and climate conditions. The commonly used variables of individual scales include tree species, tree size (a proxy for tree age measured by diameter and height, and crown width), and tree growth, such as diameter increment, basal area increment, and relative basal area increment (Monserud and Sterba 1999; Wyckoff and Clark 2000; Adame et al. 2010). Tree size, which is influenced greatly by multiple factors, could help in the search for the causes of tree death in modeling tree mortality, and was the most popular factor in previous studies (Monserud and Sterba 1999; Yang et al. 2003; Crecente-Campo et al. 2009). Moreover, stand variables related to tree mortality include stand characteristics (such as canopy cover, dominant species, mean diameter, and soil conditions), and competition indicators (Wimberly and Bare 1996). According to the characteristics and position of the subject tree and its neighbors, competition indicators are often classified as distance-independent competition indices (e.g., basal area of neighboring trees larger than the target tree) and distance-dependent competition indices (Wimberly and Bare 1996). Moreover, competition was often the important variable in predicting tree mortality in previous studies (Crecente-Campo et al. 2009; Ruiz-Benito et al. 2013). However, it still featured very high variation in predicting tree mortality due to the fact that ecosystems are complicated. In addition, climatic conditions are fundamental drivers in shaping species distribution and plant performance, and the increasing frequency of climatic extremes such as climate warming, freezing, and drought may cause physiological stress and thus lead to large-scale mortality events (McDowell et al. 2008; Zheng et al. 2012; Ruiz-Benito et al. 2013). Many previous studies have focused mainly on the variables of individual and stand scales in developing tree mortality models (Monserud and Sterba 1999; Eid and Tuhus 2001; Yang et al. 2003; Crecente-Campo et al. 2009), which may lead to limitations if the effects of climate variations are not considered.

Several studies have attempted to model the probability of tree mortality with different methods, including Weibull function (Somers et al. 1980), Richard’s function (Buford and Hafley 1985), and gamma (Kobe and Coates 1997) and neural (Hasenauer et al. 2001) networks, but the results showed little improvements over logistic regression models (Monserud and Sterba 1999; Eid and Tuhus 2001; Ruiz-Benito et al. 2013). In this study, we compared such models by classifying the independent variables into individual, stand, and climate scales with the logistic regression technique using a long-term forest inventory dataset collected from more than 2500 forest permanent plots in Sichuan Province, China. The objectives of this study were to: (1) develop individual tree models for estimating non-catastrophic tree mortality for fir, oak, and birch species; (2) select a “best” model for each of the species; and (3) evaluate how three sets of variables affect tree mortality for fir, oak, and birch.

2 Material and methods

2.1 Data

Data for this study were collected from 287 permanent plots in the Sichuan Forest Inventory (SFI) system established in 1979. The plots have been re-measured systemically every 5 years since 1988, and the most recent inventory was conducted in 2012. Each plot features a square with an area of 667 m², and the distance between plots is 4 or 8 km. Tree locations and diameters at breast height (D) greater than 5 cm were recorded. Each recorded tree was numbered and tagged for assisting subsequent inventory. All dead trees, including fallen logs, snags and stumps, were recorded in the SFI system. We used the data of plots collected in 1997, 2002 and 2007. The tree mortality in this study, therefore, refers to 5-year tree mortality.

According to the “manual of forest inventory”, all catastrophic events were categorized into four types: insects and diseases, fire, climate and geological hazard, and others. In order to model non-catastrophic tree mortality, we first excluded plots with catastrophic events. Then we excluded all the managed plots and the natural-origination plots with cutting trees exceeding 5 % of both numbers and basal area during the 5-year study period. Finally, 287 plots with 8,684 trees, including 223 dead trees, were selected for modeling non-catastrophic tree mortality. Most of the plots were located at the Tableland and alpine valleys in the Midwest Sichuan Province, with just a few plots located at the east around the edge of Sichuan Basin (Fig. 1). The climate in the study area featured concentrated rainfall in summer, and mild and cloudy weather in winter, which was highly heterogeneous and influenced strongly by the Asian monsoon climate.

Fig. 1

The spatial distribution of the 287 plots in the study area

The current study focused mainly on three key species: fir (living trees: 3,642, dead trees: 59), oak (living trees: 2,731, dead trees: 68) and birch (living trees: 2,088, dead trees: 96), and the 5-year tree mortality for fir, oak and birch was 1.59, 2.43, and 4.40 %, respectively. Fir (Abies faxoniana Rehd. et Wils.) is an evergreen coniferous species with good shade tolerance and slow growth, suited to wet and cold weather and well drained soil with a height of about 40 m, which was always mixed with birch, spruce, and oak. Oak (Quercus aquifolioides Rehd. et Wils.) is an evergreen broad leaf species with a height of about 20 m, grown on dry sunny slopes or on the top of the mountain, and always mixed with birch, fir and some other broadleaf forest trees. Birch (Betula albo-sinensis Burk.) is a deciduous broad leaf species, suited to wet, bright-light weather, and always living with fir and oak. In this region, more than 80 % plots were mixed-species stands and there were also some other species, such as Picea Dietr., Larix Mill., Cupressus L. and Populus L. The three species selected for this study are distributed widely in Sichuan Province and account for about 34 % of total forest area and 43 % of standing timber volume in the province.

2.2 Variable selection

Based on the collected variables related to tree death, we divided variables into three sets: (1) individual tree variables (I) include: relBAI and diameter at breast height (D); (2) stand variables (S) include: distance-dependent competition indices (DCI), basal area of neighboring trees larger than the target tree (BAL), canopy cover (CC), stand mean diameter (SMD), soil texture include sand and clay (SAND, CLAY), C:N ratio, and total nitrogen (TN); (3) climate variables (C) include: average annual monthly precipitation in January (Pre1) or July (Pre7); average annual temperature (Tavg); average annual maximum temperature in July (Tmax); average annual minimum temperature in January (Tmin). We also included the log-transformed and square terms of these variables in building the models considering possible nonlinear relationships. The mean values and the amplitude according to these variables for each species are given in Table 1.

Table 1 Values (mean [min, max]) of the variables for live and dead trees for fir, oak, and birch. D Breast-height diameter, relBAI relative basal area increment, BAL basal area increment, DCI distance-dependent competition indices, SMD stand mean diameter, CC canopy cover, PRE average annual precipitation, Pre1 and Pre7 average annual monthly precipitation in January or July, Tavg average annual temperature, Tmax average annual maximum temperature in July, Tmin average annual minimum temperature in January, TN, total nitrogen, CN ratio of soil organic carbon to nitrogen, SAND and CLAY fractions of soil texture

2.2.1 Individual-tree variables

Relative basal area increment (relBAI) was calculated as the ratio of BAI (basal area increment) and BA (basal area) of individual trees for the growth-dependent mortality model (Bigler et al. 2004). In contrast to other growth (Wyckoff and Clark 2000) or growth trend (Bigler and Bugmann 2003) variables, the variable of relBAI takes into account the relative change in tree size. In this study, relBAI (Wunder et al. 2007) for each tree i, was calculated by dividing the BAI (1997–2002) by its basal area (BA) in 2002:

$$ \mathrm{relBAI}\left(i,{t}_1,{t}_2\right)=\frac{\mathrm{BA}\mathrm{I}\left(i,{t}_1,{t}_2\right)}{\mathrm{BA}\left(i,{t}_2\right)} $$

(1)

In this equation, relBAI (i, t ₁, t ₂) represents the relBAI of tree i during the years t ₁ to t ₂, with its value restricted to the range of 0–1; BA (i, t ₂) is the BA of tree i in the year t ₂; BAI (i, t ₁, t ₂) is calculated by BA of tree i in year t ₂ minus BA of tree i in year t ₁; t ₁ is 1997 and t ₂ is 2002.

This reflects the relative increment of tree size; thus, small trees have a higher relBAI for the same BAI. Furthermore, we selected trees from the forest inventory dataset re-measured in 1997, 2002 and 2007, when all trees should be alive at the first and second inventory but might be dead by 2007. Since we did not know exactly the year when trees died between 2002 and 2007, the growth time for dead trees would be less than 5 years and, thus, the relBAIs (i, 2002, 2007) of trees were not used. Diameter at breast height was used in lieu of tree age (Adame et al. 2010). The nonlinear or U-shaped curve effects can be captured with a hyperbolic D⁻¹ or D and D² transformation of diameter, respectively (Monserud and Sterba 1999; Adame et al. 2010).

2.2.2 Stand variables

We chose a widely used distance-independent competition index of BAL (sum of BA of neighboring trees larger than the target tree) (Contreras et al. 2011) to fit our models. The competition index was calculated as:

$$ {\mathrm{BAL}}_i={\displaystyle \sum_{j=1}^{\mathrm{n}}}\ \uppi \cdot {\left(\frac{{\mathrm{D}}_j}{2}\right)}^2,\left({\mathrm{D}}_j\in P\right),P=\left\{{\mathrm{D}}_j,j:1,2,\dots, n\Big|{\mathrm{D}}_j>{\mathrm{D}}_i\right\} $$

(2)

Where, D _j is diameter of tree j, which is larger than that of the target tree i. BAL can explain the competition relationship between trees of different diameter classes. An increasing value of BAL leads to a reduction in diameter growth. Larger trees typically have a more competitive advantage for limited resources such as light, water and nutrients over smaller ones, which have little effect on the growth and survival of larger trees (Yang et al. 2003).

Within the square plots, two-dimensional Cartesian coordinates (x _i, y _i) were used to define each tree’s position, which can be assumed as a spatial point process ignoring tree size (Penttinen et al. 1992). Qualitative (e.g., species) and quantitative (e.g., D, height) values are often attached to these points. Thus, another distance-dependent competition index was calculated to fit our models by the position and diameter of trees. However, edge effects may have effects on spatial traits calculations, thus a buffer zone could be set around the plot giving an unbiased index, but this would also exclude the edge trees of the plots (Rouvinen and Kuuluvainen 1997). In response, a widely used translation edge-correction method was applied that copied the same pattern of plots to its adjacent eight directions (Pommerening and Stoyan 2006), and the results were then used to calculate a distance-dependent competition index (DCI). The sum of horizontal angles distance-DCI (Rouvinen and Kuuluvainen 1997) was calculated as follows:

$$ \mathrm{D}\mathrm{C}\mathrm{I}={\displaystyle \sum_{i=1}^n}\left({\mathrm{d}}_{\mathrm{i}}/{\mathrm{d}}_{\mathrm{j}}\right)\times \arctan \left(\frac{{\mathrm{d}}_{\mathrm{i}}}{{\mathrm{d}\mathrm{ist}}_{\mathrm{i}\mathrm{j}}}\right) $$

(3)

In this equation, d _j is the diameter of the subject tree, d _i is the diameter of the competitor tree; dist _ij is the distance between the subject and the competitor trees. If a subject tree has more large neighbors, it will be subject to a higher competitive stress.

Other stand scale variables, such as canopy cover (CC) and stand mean diameter (SMD) were obtained directly from the SFI dataset. Soil texture (SAND, CLAY), C:N ratio, and total nitrogen employed to explain soil aeration, available water and nutrients, were extracted from the 1:1,000,000 soil database of the Second National Soil Survey of China (Shi et al. 2002), which integrated 94,000 features and 7,292 soil profiles.

2.2.3 Climate variables

The regional climate variables including, average annual temperature (Tavg) and precipitation (Pre), average annual monthly precipitation (Pre7) and maximum temperature (Tmax) in July, average annual monthly precipitation (Pre1) and minimum temperature (Tmin) in January, were derived from China’s ground basic meteorological data products (1997–2002). We employed the widely adopted thin-plate splines method of ANUSPLIN to interpolate the weather station data to each plot with a resolution of 500 × 500 m (Hutchinson 1998). The R square values between observation and prediction through an independent validation dataset were 0.95–0.98 for temperature and 0.8–0.85 for precipitation.

2.3 Model development

Logistic regression models were developed for each of the three species to estimate the probability of survival using explanatory variables described in Table 1. The logistic regression models were fitted by maximum log-likelihood method using R (R version 3.0.2; Ihaka and Gentleman 1996). The backward stepwise selection method was applied to generate models utilizing only significant variables. It was convenient to fit the probability of survival of each individual tree with a binary response variable (living or dead, 1 or 0), which was derived from the survival probability P = P(Y = 1|x ₁, x ₂ … x _n), and can be described as:

$$ {\mathrm{P}}_{\mathrm{s}}=\frac{{\mathrm{e}}^{\upalpha +\upbeta \cdot \mathrm{X}}}{1+{\mathrm{e}}^{\upalpha +\upbeta \cdot \mathrm{X}}} $$

(4)

Thus, the “logits” sums as:

$$ \log \left(\frac{{\mathrm{p}}_{\mathrm{s}}}{1-{\mathrm{P}}_{\mathrm{s}}}\right)=\upalpha +\upbeta \cdot \mathrm{X} $$

(5)

Where, P _s is the probability of survival; X is the matrix of predictor variables used to fit the model; α is the constant parameter values; and β is the vector of coefficients. Therefore, we can fit the survival probability as a linear regression model (McCullaugh and Nelder 1989).

We built the models with combinations of the three scales of variables, such as I, S, C, I + S, I + C, S + C, I + S + C, and then selected the best one in each combination. Specifically, we first built the full models with all the variables in each combination, and then got rid of the variables when the multicollinearity occurred (the models often containing too many variables plus the transforms of variables), using Akaike information criterion (AIC) values to drop variables when too many variables entered in case of overestimated (using “drop1” method in R). Furthermore, we kept using the backward stepwise selection procession in each modeling and tried all possible combinations of the significant variables at each level. All the models were named with species plus one of the seven compositions of three scales plus the serial number of each model (1 – n), such as FirI1 – FirIn, FirIS1–FirISn, FirISC1–FirISCn, etc.

2.4 Model evaluation

The logistic model can be evaluated in many different ways. AIC is a statistic commonly used to evaluate model performance. First, we used a relative AIC (△AIC), which is the difference between the AIC of the current model and the lowest AIC of all models, where a lower △AIC value indicates a better model fitting (Bigler and Bugmann 2004; Wunder et al. 2008). We selected the “best” models with lower △AIC values for each species. Second, we applied the receiver operating characteristic (ROC) curve plotting the true positive rate (sensitivity) against the false positive rate for different cut-off values to indicate the discrimination power of different models (Fielding and Bell 1997). The area under the ROC curve (AUC) is a threshold-independent indicator to differentiate live and dead trees. An AUC value of 1.0 means all the live and dead trees are perfectly predicted, yielding a point (0, 1.0) on the y-axis in the ROC plot. The closer the ROC curve toward the upper left corner in the ROC plot the better discrimination between live and dead trees (Zweig and Campbell 1993). In general, an AUC value of 0.5 means no discrimination, 0.6 – 0.7 indicates poor discrimination, 0.7 –0.8 means acceptable discrimination, and 0.8 –0.9 means excellent discrimination between live and dead trees (Crecente-Campo et al. 2009).

Furthermore, we evaluated the models using Hosmer-Lemeshow’s (H-L) goodness-of-fit statistic method, which assessed how the observed tree mortality matched with the predicted mortality (Hosmer and Lemeshow 1989). All the sample trees for each species were sorted in an ascending order according to the predicted probability and the trees were divided evenly into ten groups based on the predicted probabilities (Yang et al. 2003). The Chi-square statistic was employed to test the difference between the observed and modeled mortalities. Here we used both Chi-square value and corresponding P-value of the H-L goodness-of-fit statistic to assist our model evaluation and selection. P-values less than 0.05 mean that the modeled tree mortality differed significantly from that of the observation. On the other hand, a greater P-value means that the modeled mortality was closer to the observed mortality.

In statistics, the variance inflation factor (VIF) quantifies the severity of multicollinearity in an ordinary least squares (OLS) regression analysis. A lower VIF value is expected with all other factors being equal. The VIF equals 1 when the current vector is orthogonal to each column of the designed matrix. However, VIF is usually larger than 1, and an accepted VIF value may be less than 10, 5 or even 4 (Pan and Jackson 2008). We adopted the accepted VIF value of 4 in this study. We used the R statistical software to perform all calculations (R version 3.0.2; Ihaka and Gentleman 1996).

2.5 Model validation

We used the commonly accepted k-fold cross-validation application to validate our models (Harrell 1998). A single subsample out of k subsamples is retained as the validation data for testing the model, and the remaining k-1 subsamples are used as fitting data. The process will repeat k times to make sure that each subsamples are used exactly once as the validation data. In this study, the original sample was randomly classified into 10 equally sized subsamples. We calculated AUC values for the validation based on the predicted survival probabilities.

3 Results

3.1 Model development and evaluation

Our results showed that FirISC1, OakISC1, and BirchISC1 were the “best” models for simulating tree mortality for species of fir, oak, and birch (Table 2). For the fir species, the individual-based model FirI1 was improved slightly by adding stand variables of DCI and CC to the model (FirIS1), and was improved significantly by adding climate variables of Tavg, Pre1 and Pre. Furthermore, the model was improved greatly by using individual, stand, and climate variables simultaneously. The “best” model, FirISC1, was achieved by including the relBAI’, DCI, CC, Tmax, Pre1, and Pre as the independent variables, the corresponding AIC value was 0.0, and AUC value increased to 0.795. The “best” model FirISC1 also had a lower Chi-square value of 5.4 and the corresponding P value was 0.71. The P value was much larger than 0.05, which suggested that the model FirISC1 was a very good fit of the data (Table 2). Similarly, we also found that the model performance was improved significantly by using the three levels of variables among all the models for oak and birch trees (Table 2). The “best” model of the two species, OakISC1 and BirchISC1, had the lowest AIC values and the highest AUC values of 0.845 and 0.821 for the fitted data among all the models. Additionally, the Chi-square values for oak and birch were 9.25 and 11.41, and the corresponding P values were 0.32 and 0.18, which means the two models were good fits of the data as well.

Table 2 Model comparison and selection for fir, oak, and birch. AIC Akaike information criterion, AUC area under receiver operating characteristic curve

We examined model multicollinearity using the VIF method and found that all VIF values ranged from 1.13 to 1.60, 1.15 to 2.23, and 1.06 to 2.33 for each variable of the “best” three models for fir, oak, and birch, respectively. This means that multicollinearity is not a big concern in our “best” models. Moreover, VIF values for variables of all the models for fir, oak, and birch ranged from 1.00 to 2.03, 1.00 to 3.62, and 1.01 to 2.95, respectively, except that the square term of D and D², Pre1 and Pre1², and Tmax and Tmax² with VIF was larger than 10. The correlation matrix indicated that the correlation coefficients among climate variables was much higher than that of other variables. For the individual and stand scales of variables, most correlation coefficients were lower than 0.5 (Supplementary Table S2a, S2b, S2c).

3.2 Model validation

We validated all models for the three species using a 10-fold cross-validation method. The mean validated AUCs for fir, oak, and birch were 0.690 (range 0.510–0.780), 0.727 (range 0.560–0.828), and 0.725 (range 0.554–0.810), respectively, and the mean fitted AUC for fir, oak, and birch were 0.716 (range 0.555–0.795), 0.733 (range 0.574–0.845), and 0.741 (range 0.583–0.821), respectively. Moreover, the validated AUC values of cross-validation data for all models of fir, oak and birch declined slightly by 3.63 %, 0.82 % and 2.16 %, respectively. In addition, the validated AUC values for the “best” models of fir, oak, and birch were 0.780, 0.828, and 0.810, respectively, which were slightly lower than the AUC values of fitted data, but still were the highest values among the validated values of all the models.

3.3 Factors affecting tree mortality

The significant variables for all the three species included relative basal area increment (relBAI), Diameter at breast height (DBH), distance-dependent competition index (DCI), canopy cover (CC), average annual maximum temperature in July (Tmax), and average annual monthly precipitation in January (Pre1).

We found logarithm-shaped relationships between relBAI and survival probabilities for fir, oak and birch with models built with a single variable of log-transformed relBAI (relBAI’). The coefficients of relBAI’ were positive for fir, oak and birch, indicating that the probability of tree mortality may be less when the relBAI is large. The mean values of relBAI of fir, oak and birch for the living trees were 0.144, 0.136 and 0.170, respectively, decreasing to 0.117, 0.101 and 0.105, respectively, for dead trees (Fig. 2a). Moreover, for a slow-growing tree, a small increase in growth led to a relatively large decrease in death probability, while for a fast-growing tree, a large increase in growth led to a relatively small decrease in mortality for fir, oak and birch (Fig. 2b). In addition, diameter was also a significant predictor for the three species. The positive coefficient for diameter (D) and the negative coefficient for diameter squared suggested a U-shaped mortality trend for fir and oak (Table 2). However, the coefficients of diameter (D) were positive for birch without the squared term of diameter, indicating that tree mortality would be less as diameter increased.

Fig. 2

a Relative basal area increment (relBAI) for fir, oak, and birch measured between 1997 and 2002. Box plots The box length covers the first and third quartile value of the trees, and the bold line inside the box is the median. Notches define a confidence interval of the median. Whiskers extend to the most extreme data point between the box and 1.5 times the box length. b Growth–mortality (expressed in logits) relationship according to the variable relBAI for fir, oak and birch

The coefficients were negative for competition indices (DCI and BAL) for fir, indicating that the probability of tree mortality would be higher when the competition is higher. We also found that the smallest trees (less than 10 cm) endured higher competition than larger trees for all the three species in the current study (Fig. 3). There were other stand variables that significantly affect tree mortality including CC squared (Fig. 5) and SMD (or SMD squared) with positive coefficients for oak and birch, TN and SAND with positive coefficients for oak, and CN and CLAY with negative coefficients for birch. In addition, the positive TN and negative CN suggested that the tree mortality would be lower as the total nitrogen increased.

Fig. 3

Distance dependent competition index (DCI) and distance independent competition index (BAL) for fir, oak and birch distributed by diameter class

Our results also showed that tree mortality increased with the increases in Tmax (average annual maximum temperature in July) and Pre1 (precipitation in January) for fir (models FirC3 and FirC4, in Table S1a) and oak (models OakC2 and OakC3, in Table S1b), while mortality increased for birch when Tmax was higher than 21.3 °C and Pre1 was larger than 6.9 mm (by models of BirchC3 and BirchC4; Table S1c, Fig. 4). Moreover, we found that all the three species showed U-shaped responses to temperature and precipitation and these responses can be well quantified through quadratic equations (Fig. 4, R ² = 0.95 – 0.99).

Fig. 4

Relationship between tree mortality and climate variables such as a precipitation in January (Pre1) and b average annual maximum temperature in July (Tmax) for fir, oak, and birch. Tree mortality increased with the Tmax and Pre1 rising for fir and oak, respectively. Moreover, a U-shaped relationship was found for birch of Tmax and Pre1. The inflection values of Pre1 and Tmax for birch were 6.9 mm and 21.3 °C, respectively

4 Discussion

4.1 Improvement of tree mortality models

Modeling individual tree mortality is tremendously challenging due to the multiple factors and processes involved in tree senescence and the usually small proportion of dead trees in forests (Ruiz-Benito et al. 2013). Compared to models built with individual and stand levels of variables, the improved models built with both biotic and abiotic factors had better performances in predicting tree mortality. Our results showed that the improved models featured AUC values of 0.80, 0.85, and 0.82 for fir (FirISC1), oak (OakISC1), and birch (BirchISC1), respectively, which were similar to AUC values of 0.62 to 0.87 found in the most recent models reported (Wunder et al. 2008; Crecente-Campo et al. 2009; Adame et al. 2010).

Our results showed that VIF values of all variables for the three species were lower than 4.0, suggesting that the improved models did not have strong multicollinearities (Pan and Jackson 2008). Moreover, degradation in performance due to over-fitting can be surprisingly large, which might result from a higher number of predictors, even with excellent discriminatory power. In tree mortality modeling, it seems prudent to keep the number of explanatory variables below one-tenth of the number of dead trees for each species (Wunder et al. 2007). In this study, the number of predictors in the best models for all three species did not exceed this recommended ratio. Furthermore, cross-validation provides a simple and effective method for both model selection and performance evaluation, which can help reduce over-fitting during model selection (Cawley and Talbot 2010). The AUC values of the cross-validation and fitted data were very close for the “best” models, suggesting an unbiased validation of the models. Therefore, the improved performance of our models was unlikely to have arisen from the over-fitting of the models.

Many previous models were developed using variables at individual scale, such as basal area and tree diameter (Monserud and Sterba 1999; Bigler and Bugmann 2004; Wunder et al. 2007). However, an earlier study found that individual level variables combined with stand age could improve tree mortality modeling in jack pine forests in Kimberly Clark in Canada (Chen et al. 2008). Moreover, a recent study found that tree mortality was affected strongly by competition, annual precipitation and temperature (Ruiz-Benito et al. 2013). Interestingly, in the current study we also found that adding climate scale of variables could further improve the models and provide explanations of abiotic relationships with tree mortality. However, it should be noted that climate elements are probably highly correlated with geographic variables such as latitude, longitude, and elevation, which may confound climate factors (Xu et al. 2004). Moreover, we plotted the residuals versus plot locations and elevations and found that the residuals might have some spatial structures (Fig. S1). In addition, some plots had extremely high residuals, which was due mainly to very few trees in some plots (for example, some plots had only one tree). However, we also found that the residuals had no significant relationships with latitude, longitude, and elevation (Table S3). Therefore, in this study we excluded geographic variables and instead employed variables of climate scales for developing the models.

4.2 Effects on tree mortality of individual and stand factors

Previous studies have found that the relBAI of dead trees was considerably lower than that of surviving trees in the same stand in forest inventory datasets (Wunder et al. 2007, 2008; Castagneri et al. 2010) and tree ring data (Bigler et al. 2004). In this study, we also found that the growth rate of surviving trees is much higher than that of dead trees for all three species (Fig. 2a), and the log-transformed relBAI made a significant contribution to tree mortality modeling (Fig.2b, Table 2). Earlier studies also found that the probability of survival for different tree species increased with increasing log-transformed relBAI, as hypothesized by Waring and Pitman (1985). Such a correlation was also found in seven species of unmanaged forests in European spruce (Wunder et al. 2008) and Norway spruce in the eastern Italian Alps (Castagneri et al. 2010).

Usually, the total photosynthetic surface area of a stand expands as the stand develops and reaches a maximum level and then stabilizes or declines due to nutrient limitation, while nonphotosynthetic tissues increase continuously (Smith and Long 2001). As a result, the ratio of photosynthetic to nonphotosynthetic tissue declines as trees grow larger, which may result in increased maintenance costs and subsequently less photosynthate available for tree growth (Yang et al. 2003). In this study, we found that the mortality of fir and oak shared a U-shaped mortality curve, with an increased mortality rate for large old trees (Model FirI2 in Table S1a and OakIC1 in Table 2), suggesting that these large old trees lacked vigor and suffered lower survival with increasing diameter. Similar results were also found in Norway spruce (Picea abies; Monserud and Sterba 1999) and White spruce (Picea glauca; Yang et al. 2003). It has been noted that the net photosynthetic rate per unit leaf area declined for large old trees, which might reduce growth efficiency, and increase the rate of tree mortality (Waring and Schlesinger 1985). Furthermore, we found that larger trees (> 80 cm for fir, and > 50 cm for oak) were located in areas of more suitable temperature and precipitation (most were at the lower range of Tmax7 and Pre1 in Fig. 4), suggested that the higher tree mortality for large trees was more prone to relate to tree senescence.

Competition can be explained as a negative reciprocity between trees through direct interference or indirect exploitation of shared resources (Connell 1990). Trees may compete directly only with their immediate neighbors, and there is a strong correlation between neighbor distances and tree size in inter-tree competition (Getzin et al. 2006). Intensified competition between trees reduces the growth of individual trees and may lead finally to the death of trees. In this study, our results showed that dead trees experienced higher competition intensities than living trees (Table 1). Interestingly, we also found that competition indices were important for improving mortality models. A study by Rozas and Garcia-Gonzalez (2012) found a high correlation between oak tree mortality and inter-tree competition using a distance DCI.

Soil texture is highly related to soil moisture and soil fertility, which are critical to tree growth, and the resistance to various stresses may affect tree mortality, especially for light-demanding and fast-growing pioneer species (de Toledo et al. 2011). In our study, we found that tree mortality decreased with a sandier soil texture for oak and lower clay content for birch. Moreover, higher sand and lower clay content in soil related to higher total nitrogen content (Table S2b, S2c), which may suggest that well-drained fertile soil benefits oak and birch trees. Similarly, Fensham and Fairfax (2007) found that tree death was greater as the clay content increased. However, Davies (2001) found that trees had a lower growth rate when the soil texture turned sandier in the mixed dipterocarp forest in Malaysia.

The forest canopy structure affects the light distribution, influences local stand precipitation, air humidity and movements, and even soil moisture (Jennings et al. 1999). The canopy affects plant growth and survival, hence determining the nature of the vegetation and wildlife habitat. The results of this study show that tree mortalities decreased as the canopy cover increased for fir and oak, suggesting that increased canopy coverage favors more light for dominant and constructive species (Fig. 5). However, we also found that birch shared an inverse-U-shaped curve, suggested that shade and light are not the dominant factors for the pioneer species of birch in a closed canopy, and the lower tree mortality at the beginning stage of succession might be due to the availability of sufficient light under an open canopy. Moreover, the tree mortality of birch was the highest when the stand canopy cover was about 0.65 due to the large number of dead small trees (DBH < 10 cm).

Fig. 5

Relationships between tree mortality and canopy cover for fir, oak and birch

4.3 Effects of climate conditions on tree mortality of three species

The individual species studied here thrived only within definite ranges of environmental conditions according to ecological niche theory (Hirzel and Le Lay 2008). The most suitable niches of a species are located in the middle of its distribution range where the temperature and precipitation are moderate. In agreement with this theory, our results suggested that the species with most climatically suitable conditions had the lowest tree mortalities. The high mortality at both extremes (e.g., the higher and lower ends of Tmax) was attributed to various biotic and abiotic interactions, such as more frequent freezing, drought, insect and disease damage, and more severe competition with other species in the marginal areas. In our case, birch trees had better optimal precipitation (January) and temperature (Tmax) than fir trees, which had better optimal values than oak trees (Fig. 4). Moreover, birch had the highest tree mortality partly because it had a narrower ecological amplitude and was a pioneer species in a community that colonized rapidly at the beginning of the succession stage with a short life history. However, as a dominant species, fir was shade-tolerant and had a long life history, had a lower tree mortality than birch and the middle succession species of oak, which had the most moderate tree mortality of the three species. It was noted that the oak species in the current study is an alpine oak, which is distributed in very dry and cold alpine regions.

Climate-driven tree mortality can be affected by both heat stress and increases in frequency, duration, and/or severity of drought (Williams et al. 2013; Zhang et al. 2014). The heat-induced stress and drought under climate change could further change the composition, structure, and biogeography of forests in many regions (Allen et al. 2010). Moreover, climate-induced stress is associated with other processes such as insect outbreaks and wildfire, which may potentially increase tree mortality (Brando et al. 2014). In this study, we found that maximum July temperature significantly raised tree mortality of the three species during 1997–2002, which included one of the warmest years (1998) in recent times in Southwest China (Wang and Gong 2000). Many previous studies also found that the rates of tree mortality increases with elevated temperature (Van Mantgem et al. 2009; Williams et al. 2010; Ruiz-Benito et al. 2013; Zhang et al. 2014). Furthermore, increased summer temperatures or rapid climate warming have had a strong negative effect on tree growth (Dulamsuren et al. 2013; Liu et al. 2013), and thus have led to increased tree mortality in many regions (Williams et al. 2010).

In addition to temperature, climate warming-induced drought may be another dominant driver of the widespread increase in tree mortality (McDowell et al. 2008; Williams et al. 2013). Several previous studies have found that drought was one of the dominant drivers for the increase in tree mortality (Van Mantgem et al. 2009; Peng et al. 2011). Our results suggested that the predicted warming temperature and declining precipitation across large areas in Sichuan Province until 2070 (Hijmans et al. 2005) might potentially exacerbate severe tree mortality, thus leading to the loss of carbon sequestration for the widespread species of fir and oak that are the main carbon pools of forests.

Accumulating snow on tree crowns and stems in forests may cause freezing damage, stem and branch breakage, even causing trees to fall down, which is considered as an important factor related to tree mortality (Nykanen et al. 1997; Zheng et al. 2012). Winter precipitation is mostly snow or freezing rain in the studied area, which may affect the performance of trees, especially at high elevations (Walsh et al. 1994). The characteristics of taper and crown are the main factors controlling the resistance of trees to snow (Walsh et al. 1994). Future winter precipitation is predicted to increase in Sichuan Province (Hijmans et al. 2005), which may lead to increased tree mortality for fir, oak and birch. However, the predicted rise in winter temperature might help to mitigate the effects of freezing. Moreover, the prolonged duration of snow cover in spring reduces the growth season, which may also be crucial for the survival and successful growth of trees in the subalpine timberline (Senn 1999).

5 Conclusions

We developed logistic regression models to predict individual-tree mortality for three tree species: fir, oak, and birch. Our results showed that including climatic variables improved non-catastrophic tree mortality models for the three species. The selected models featured a good level of discrimination, suggesting that tree mortality can be better predicted by combining the variables of individual, stand, and climate scales together in Southwest China. Moreover, the models are also suitable for predicting stand dynamics and thus forest ecosystem structure and function in Southwest China due to the fundamental components of individual-based growth and yield models. In addition, future climate change such as climate warming, drought, and freezing events may threaten the survival of these three tree species.