© INRA, EDP Sciences, 2010

1. Introduction

Forest growth models have become an indispensable tool for forest management (Vanclay, 2003). Tree mortality is an important process in forest ecosystem dynamics and is one of the most complex phenomena in forest modelling. Severe inter-tree competition throughout stand development gradually decreases the growth of individual trees and finally leads to the death of some trees (Jutras et al., 2003). This process tends to be gradual but can occasionally be abrupt (Waring, 1987). Modelling this mortality is difficult because of the complex interactions between different environmental stresses, minimal understanding of whole-plant mortality processes, and a chronic shortage of data (Hawkes, 2000). However, it is fundamental to include a mortality function in forest growth models since tree mortality produces an open space which leads to (Franklin et al., 1987): (1) increased growth of surrounding trees through reduced competition, (2) opportunity for regeneration by creation of gaps, and (3) the addition of large woody debris to the structure of stands.

Natural mortality can be separated into regular and irregular mortality. Regular mortality, or self-thinning, is due to competition for light, water and soil nutrients within a stand (Peet and Christensen, 1987). Irregular mortality results from random disturbances or hazards such as fire, wind, snow or insect outbreaks, where the probability of survival is very low at the tree and stand level but not necessarily at the landscape level (Alenius et al., 2003). Generally, when modelling tree mortality or survival probability, only regular mortality is considered (Amateis et al., 1997; Monserud and Sterba, 1999).

A variety of empirical approaches have been used for describing and modelling tree mortality at different scales. Stand level mortality models predict stand density changes over time and they often estimate mortality in terms of volume, basal area or number of trees per hectare, while individual tree mortality models predict the probability of survival or death for each individual tree involved in growth projections. Stand averages in terms of mortality will be less useful because of the high degree of variability (Pedersen, 2007), but the scale selection also depends on the other components of the growth model. Individual tree models, though more complicated, are often preferred from a management perspective (Vanclay, 1994), because they provide a more detailed description of the stand structure and its dynamics (Mabvurira and Miina, 2002).

Figure 1 Map of sample locations and biogeoclimatic strata (Elena Roselló, 1997).

Several statistical methods have been used to develop empirical individual tree mortality models: logistic function (Monserud and Sterba, 1999), two (Álvarez González et al., 2004; Diéguez-Aranda et al., 2005) and three (Fridman and Stahl, 2001; Meng et al., 2003) step approaches and neural networks (Hasenauer et al., 2001; Weingartner et al., 2000). In individual-based models, tree death is predicted based on a range of variables that can be classified roughly as (Bigler and Bugmann, 2003; Eid and Tuhus, 2001; Monserud and Sterba, 1999): (i) size-related variables such as diameter at breast height or tree height, (ii) growth-related variables such as competition and measures of stand density, (iii) measures of tree vitality.

Data from National Forest Inventories (NFI) with permanent plots probably provide the best data set to reflect the full range of variability with respect to treatments, sites, forest structures, and tree species. The disadvantages of using this type of data are first, that the methodology is not specifically designed with model development, and second, the difficulty involved in identifying the dead trees. Models for individual tree mortality based on such data have previously been developed for different species and regions (Eid and Tuhus, 2001; Hasenauer et al., 2001; Monserud and Sterba, 1999). The Spanish NFI has been used for modelling tree mortality in pine species (Trasobares et al., 2004a; 2004b).

The present analysis focuses on pyrenean oak (Quercus pyrenaica Willd.), a widespread tree species in the northwestern part of the Iberian Peninsula. The increasing interest in the use of this species for either direct production (such as wine barrels or biomass) or indirect production (such as silvopastoral uses, recreation or environmental preservation) justifies the need to develop a dynamic growth model (Adame et al., 2008a; 2008b; Cañellas et al., 2004).

A number of individual tree mortality models have been developed for oaks growing in Central and Northern Europe (Eid and Tuhus, 2001; Fridman and Stahl, 2001; Monserud and Sterba, 1999). All of them used data from National Forest Inventories and species of the genus Quercus were considered to be a homogeneous group or even included within “other deciduous” species. Red and white oak groups in the eastern part of North America have also been studied (Shifley et al., 2006). However, a mortality model for Mediterranean oak, like Q. pyrenaica has, to date, not been developed.

The purpose of this paper is to model 10-year individual mortality in pyrenean oak for a large data set belonging to the Spanish National Forest Inventory. Only the regular mortality caused by competition, old age, and small-scale attacks from insects and diseases were modelled in this study. The main objectives are to identify which site/ stand characteristics or tree variables are the principal factors associated with oak mortality, to assess the strength of the relationships and interactions among variables and to develop a mortality function to complete the individual tree model for pyrenean oak.

2. Material and methods

2.1. Data set

Data for Pyrenean oak stands from the Second and Third Spanish National Forest Inventory (SNFI) located in the northwest Spain (Castilla y León region) were used for the purposes of this study. The aforementioned region accounts for more than 53% of the species distribution area in Spain (a total of 1208768 ha) (DGCN, 2005) (Fig. 1). The SNFI is a systematic sample of permanent plots of variable radii distributed on a 1 km square grid. The Second Spanish National Forest Inventory was conducted between 1993 and 1994 and measurements were repeated for the first time after 10 y, i.e. between 2003 and 2004 (Third SNFI). The plots, which are circular in shape, are composed of four sub-plots with radii of 5, 10, 15 and 25 m, and a minimum diameter at breast height threshold of 75, 125, 225 and 425 mm, respectively. Tree species, diameter at breast height (dbh) and height (h) were recorded for all inventoried trees. Tree mortality is a discrete event, the only possible values being 0 (alive) or 1 (dead), and was recorded by counting trees that had died between the second and third inventory. The values for basal area (BA), quadratic mean diameter (DG) and number of trees (N) were calculated by expanding individual tree measurements (Bravo et al., 2002).

Because of the high variability in silvicultural and ecological conditions in pyrenean oak stands, as well as the lack of past management information in the SNFI, a second selection was made according to a stand typology classification for Q. pyrenaica stands (Roig et al., 2007). The selection of SNFI plots excluded impoverished or very sparsely populated stands (Reineke’s stand density index under 200). Since no information was available with respect to why and how harvesting was carried out, plots that had undergone any harvesting treatment between the two measurements (according to the SNFI) were also excluded. A total of 219 permanent plots and 4505 trees were selected to fit the mortality model. At remeasurement, all trees were alive in 163 sample plots (74%) and there were 160 dead trees, corresponding to an overall annual mortality rate of 0.35% (Tabs. I and II).

Six strata according to biogeoclimatic characteristics (precipitation, temperature, altitude and soil characteristics) are identified in the studied region (Elena Roselló, 1997) which can be used as explanatory variables in the mortality model (Fig. 1). In addition, ecological parameters calculated according to the location of each plot (altitude, geographic position and catchment area) using the climatic estimation models proposed by Sánchez et al. (1999) were also evaluated (Tab. III).

Independent data sets from SNFI permanent sample plots located in another region (Madrid region) with similar biogeoclimatic characteristics were used to test the performance of the proposed model. Measurements were taken and data collected at these sample plots using the same procedure used for the modelling data set. Therefore, 504 trees from 29 plots were considered (Tab. IV).

2.2. Model development

2.2.1. Selection of variables

Selection of appropriate predictor variables should not only be based on test statistics, but also on a basic understanding of how forest ecosystems function and how factors contributing to mortality are expressed (Pedersen, 2007). Site index and age were excluded intentionally because it is not possible to establish the age of every tree.

Tree size has been used as a proxy of tree age (Yao et al., 2001). Individual diameter (dbh) and height (h) are available to characterize tree size. For most species, the mortality rate for the smallest trees is quite high, and declines rapidly as individuals survive to larger diameter classes. This nonlinear effect can be captured with a hyperbolic dbh^-1 transformation of diameter (Monserud and Sterba, 1999). The oldest individuals show signs of decline or senescence, so the probability of mortality increases, forming a U-shaped curve which could be modelled by including terms for dbh and dbh² in the exponent of the logistic regression (Monserud and Sterba, 1999).

A major mortality agent is inter-tree competition (Yang et al., 2003), which can be either one-sided or two-sided. In one-sided competition, larger trees are at competitive advantage over smaller trees and neighbours do not affect the growth and survival of larger trees (Cannel et al., 1984). One-sided competition has been expressed here by the ratio of the diameter of the subject tree to the quadratic mean diameter of the sample plot (dbh/DG) (Hamilton, 1986), the summarised basal area for all trees greater than the subject tree (BAL) (Monserud and Sterba, 1999), the ratio of BAL and basal area of the stand (BAL/BA) (Jutras et al., 2003), and the ratio of the height of subject tree to the dominant height of the sample plot (h/Ho) (Avila and Burkhart, 1992; Zhang et al., 1997).

In two-sided competition, all trees impose some competition on their neighbours, regardless of their size (Cannel et al., 1984), and this can be described by stand variables such as density and basal area. The following stand variables were tested for inclusion in the model: stand density (N, stems/ha), quadratic mean diameter (DG, cm), basal area (BA, m²/ha), dominant height (Ho, m), Hart index (HSI) and Reineke index (SDI).

Earlier studies have found that sites with different levels of productivity exhibit differences in mortality (Eid and Tuhus, 2001; Yao et al., 2001). Although site index data was not available, it was possible to include different site specific variables (Tab. III).

Independent variables to be included in the function were selected variable by variable through visual analysis of data plots, as well as through backwards logistic regressions (SAS/STAT, 2000). The set of variables should be a combination of the different groups (tree size, competition and site specific variables), avoiding correlations between them.

2.2.2. Multilevel model with binary response

The presence or absence of mortality (M) during a 10-year period, observed for a single tree, is a discrete event. i.e. M is a “0, 1” variable. The following standard logistic function was selected for the present work (McCullagh and Nelder, 1989): (1)where π is the probability of mortality (survival is given by 1 − π) binomially distributed, α is an intercept and βX is a linear combination of parameters β and explicatory variables X.

The SNFI sample plots formed a two-level hierarchy with several trees measured in each plot, so it can be assumed that the binary responses of sample plots (occurrence of tree mortality) are correlated within these levels. A generalized linear mixed approach is required to correctly address the hierarchical structure of the data and ensure unbiased testing of all model parameters. The observed tree measurements are derived from unequal probability sampling (concentric circle depending on dbh), so they were weighted.

The model specification was made according to the structure laid down by Jutras et al. (2003). The distribution for the response is y_ij  ∼  Bin , where y_ij is the observed response of tree i in plot j and π_ij is the corresponding response probability. The logit link function is obtained as: (2) where X_ij is the transposed vector of the independent fixed variables, β is the vector of the estimated parameters and u_j represents the random deviation for j plot.

Previous research has identified certain problems with the penalized quasi-likelihood (PQL) method if the plot size is small (McCulloch and Searle, 2001), so fixed effects and random parameters at plot level (u_j) were estimated using the maximum likelihood (ML) method (Alenius et al., 2003). A multilevel model was built using the SAS/GLIMMIX macro (SAS/STAT, 2001).

2.2.3. Evaluation

A two-step procedure was used to evaluate the model, which included qualitative as well as quantitative examinations. The first step was to evaluate the model fitting statistics. Statistically, the model evaluation was based on the significance of the parameter values of the candidate explanatory variables. If a variable explains the data well, the p-value associated with that variable should be smaller than 0.05, indicating significance at the 95% confidence level. The value of the stand-level variance is also analysed as a rough indication of model improvement (Jutras et al., 2003). The Receiver Operating Characteristic (ROC) curve and the area under the ROC curve (AUC) were calculated.

In step two, the characterisation of the model error was analysed, based on the independent data set testing. The model was applied to the test data set to predict the probability of mortality for each tree. The predicted and observed mortality were then compared by visually studying deviations over the explicatory variables included in the model. Three cut-points were tested; the first was the overall survival rate found for the species, the second was the point where the sensitivity-line crosses the line of the true-negative rates and the third was a random number.

3. Results

The most important explanatory variables were the inverse of tree diameter at breast height (dbh^-1) and the ratio of the height of subject tree to the dominant height of the sample plot (h/Ho). The Pyrenean oak mortality model indicates that the higher the diameter of a tree, the lower the probability of mortality. Ratio of the height also displays an indirect relationship with mortality (the mortality decreases with an increase in the values).

Evaluation

The p-value associated with all variables in the multilevel model 2 (Tab. V) are smaller than 0.05, but the p-value of the intercept is smaller than 0.1. The value of the stand-level variance was found to be lower if the explanatory variables are included in the model. Figure 2 is the ROC curve based on the sensitivity and the specificity of the test. The area under the ROC curve (AUC) was 0.725.

Figure 2 Receiver operating characteristic (ROC) curve for the mortality model (a) and sensitivity-specificity diagram (b).

An evaluation of the model using an independent data set was carried out by studying the performance of the function, calculating the probability of mortality for each tree belonging to the test data set. Figure 3 shows the predicted probability of mortality (model data set on the left side, test data set on the right side) plotted over dbh and ratio h/Ho for observed data (dead and living trees). The behaviour of the model for both, the fitting data set and test data set, is realistic and the ascribed probabilities of mortality are considerably different between observed dead trees and observed living trees.

Figure 3 Predicted 10-year probability of individual-tree mortality of pyrenean oak as a function of tree dbh and h/Ho ratio for observed dead and live trees for model data set (a) and (c) and test data set (b) and (d).

A threshold can be used to assign mortality. If the estimated probability of mortality exceeds the threshold then the tree is considered dead. The best and the most logical choice for a threshold is the average observed mortality rate for that species (Monserud and Sterba, 1999). After ordering the predicted mortality rates, the cut-off value where the predicted mortality reaches 3% of the tree population is 0.1118. Interception point of sensitivity and specificity (0.03) is used as second threshold (Fig. 2). Random numbers can also be considered (Bravo-Oviedo et al., 2006), by running the random calculation 10 times, the threshold for each observation being the average classification rates. Table VI presents the overall predictions and the best prediction of mortality for both the fitting data set and the test data set obtained using the average observed mortality rate, intersection point of sensitivity and specificity and random numbers. The prediction of the number of dead trees is always overestimated when using the intersection point of sensitivity and specificity, although it gives the most accurate percentage (65.6%), while the average observed mortality rate and random number threshold result in no overprediction, but a much lower accuracy of predicting the correct trees (22.5% and 8.7% respectively). As can be seen from Figure 4 the mortality model behaves well with either the intersection point of sensitivity and specificity or random number, matching the observed mortality rates reasonably well for both the fitting data set and the test data set. Mortality model behaviour with average observed mortality rate seems to be slightly worst, particularly versus h/Ho ratio.

Figure 4 Relative percentages of observed and predicted 10-year mortality (according to average observed mortality rate, intersection point of sensitivity and specificity, and random number) vs. diameter at breast height (a) and (b) and h/Ho ratio (c) and (d) for model data set (left figures) and test data set (right figures).

4. Discussion

In this study, which focuses on Quercus pyrenaica, a multilevel logistic regression model was developed for predicting the probability of individual trees dead in a ten-year growth period.

The mortality process is characterised as a stochastic, rare, and irregular phenomenon (Eid and Tuhus, 2001). Stochastic individual-tree models allow the preservation of observed spatial, temporal, or nested dependencies in predictions (Fox et al., 2001). In this way, it is possible to estimate the value of the random parameters for a location not present in the original estimation data. This approach is known as localization or calibration and can be applied if supplementary observations of the dependent variable are available.

A representative sample must reflect the full range of site variability as well as the diversity of management treatments in a given population (Hamilton, 1986). This model is based on permanent sample plots measured twice by the Spanish National Forest Inventory. The data set used provides an outstanding database in terms of size (219 plots and 4505 trees) and forest conditions. This was essential since the model will be included in a growth simulator developed to cover all conditions. However, it should be borne in mind that the sampling methodology was not specifically designed to develop growth and yield models. Although the period of time (1994–2004) is considered sufficient to represent the natural mortality process of pyrenean oak, the mortality rate may have been affected by irregular climatic conditions during this period. The rate might have varied if the phenomenon had been studied over a longer period (Laard and Akça, 1997). Furthermore, the most questionable part of the sampling method concerns the determination of time of death. Although logged dead trees are recorded in the SNFI, a major source of bias is the fact that many trees may have been dead at the time of logging, but were not observed by the field crew (Fridman and Stahl, 2001).

Many mortality models for individual trees include three or fewer explanatory variables (Meng et al., 2003; Monserud and Sterba, 1999; Schröder et al., 2007). If a model is ecologically illogical, it will not perform well for any data set other than one used for model development (Hamilton, 1986), therefore all estimated coefficients were consistent with ecological expectations. The following predictor variables were included in the multilevel mortality function: individual diameter (dbh) and ratio between individual height and dominant height (h/Ho ratio). The most significant variable is the inverse of diameter at breast height, where the high mortality for small diameter trees suggests that suppressed trees are more likely to be eliminated by stand level competition. The coefficient of h/Ho is also negative, resulting in decreasing mortality rates as the values increases.

Monserud and Sterba (1999) hypothesized that the mortality rate is U-shaped, i.e. mortality decreases with increasing dbh, but increases again above a certain limit. The present model does not capture the U-shape, probably due to the low frequency of trees larger than 50 cm dbh in the data set (160 trees of 4505 total, 3.5%). This stand structure is due to the management practices applied up to the middle of the last century in this kind of stands, which were managed as either coppices or coppices with standards, by cutting trees every 20 y with subsequent shoot or root regeneration (Cañellas et al., 2004).

Crown related variables are good indicators of tree vigor (Dobbertin and Biging, 1998, 2001) and are preferred whenever possible, although in this case, such data are not available. Avila and Burkhart (1992) suggested individual height-dominant height ratio appeared to be the best variable to replace crown class variable. Mailly et al. (2009) found that depending on the species, top height trees have lower mortality rates, higher mortality rates or equal mortality rates than codominant trees or codominant and dominant trees combined, and with shade tolerant species, it appears that using codominants or a combination of codominant and dominant trees results in significantly lower mortality rates. In this study, mortality rate will be lower as individual height-dominant height ratio increases and a tree’s competitive status is more favourable.

Site variables have frequently been included in mortality models. Vanclay (1994) stated that the effect of site productivity on mortality is ambiguous: density dependent mortality should become apparent earlier in better sites, which thus display higher mortality. Better sites, however, are also capable of supporting higher densities. Moreover, several studies have shown that higher mortality is often associated with higher productivity (Eid and Tuhus, 2001; Yao et al., 2001; Zhao et al., 2007). Any bioclimatic characteristic considered in this study was significant for inclusion. This phenomenon would appear less stochastic if relevant environmental variables such as light, soil water and soil nutrients were measured on permanent plots, as well as physical aspects of microsite establishment (Davies, 2001; MacFarlane et al., 2002; Monserud and Sterba, 1999). Characterization of the genetic status of the trees would also be advantageous.

Under normal circumstances, trees will generally survive a growth period, so the percentage of correctly classified living trees is always high (Monserud and Sterba, 1999), and high prediction accuracies for dead trees are often achieved by over-predicting dead trees while reducing prediction accuracies for living trees (Crow and Hicks, 1990). This is corroborated by the average observed mortality rate and random number threshold (97.1–96.6% for fitting data set and 96.3–96.8% for test data set) and to a lesser degree, by the intersection point of sensitivity and specificity (69.3 and 58.3%). In the latter case, the percentage of correctly classified dead trees was higher (65.6 and 61.1%), although the number of dead trees is overestimated and is almost 10 times higher than observed. Dobbertin and Biging (1998) achieve prediction accuracies of dead trees between 28–36% for ponderosa pine and between 11–17% for white fir for 5-year mortality rates of 2.7 and 4.8% respectively. Dobbertin and Brang (2001) fitted a logistic regression model correctly predicted between 33–57% of the dead trees. In our study, the prediction accuracies of 22.5%, 11.1% for test data set, achieved using the average observed mortality rate (threshold of 11.18%) seems to be reasonably high given the low 10-year mortality rate in this species (3%).

Most of the published individual mortality models for oaks include the same or similar variables to those proposed in this paper. The inverse of diameter at breast height is the most important variable in the individual tree model developed by Monserud and Sterba (1999) for oak (Quercus spp.) in Austria. Eid and Tuhus (2001) determined that the probability of survival increased as diameter increased and the number of trees larger than the subject tree (BAL) decreased. Fridman and Ståhl (2001) included diameter, basal area of larger trees, altitude and stand basal area as explanatory variables for southern deciduous oaks. Shifley et al. (2006) used crown class, basal area of larger trees and individual diameter as predictors of future mortality for trees in the red oak group (species similar to pyrenean oak). Drobyshev et al. (2008) suggests that for oak, density-regulated allogenic mortality may be of higher importance than exogenous mortality, especially in the beginning of its lifespan. Drobyshev et al. (2009) found a significant and consistent association between tree condition and the probability of mortality, and significantly higher mortality rates in smaller trees of oak. The same variables as those finally included in the mortality model for Quercus pyrenaica have also been identified for others species: inverse of diameter at breast height (Bravo-Oviedo et al., 2006; Jutras et al., 2003; Karlsson and Norell, 2005); and individual height-dominant height ratio (Avila and Burkhart, 1992).

Plots that had been subjected to any harvesting operation between measurements were excluded from the data set because of insufficient information regarding treatments. Smith and Nichols (2005) compared changes in the mortality rate under a range of harvesting intensities and silvicultural systems in subtropical rainforest. Silvicultural practices, such as thinning operations, can alter mortality patterns; therefore, the mortality process needs to be monitored in stands managed with different silvicultural systems (Diéguez-Aranda et al., 2005).

In the case of pyrenean oak stands, the anthropogenic impact from villages located nearby takes the form of recurrent fires, livestock grazing, collection of firewood or poles as well as harvesting of non-timber products (Cañellas et al., 2004). In addition, the regeneration of pyrenean oak generally takes place through stump and/or root sprouting; so stand dynamics differ from those of others species. Further research into stand evolution and management alternatives is necessary given the progressive abandonment of traditional uses. The capability to associate individual tree traits with mortality has potential benefits for the study of natural processes (Zens and Peart, 2003). The improvement in tree mortality predictions allows understanding and more accurately projecting forest development.

Acknowledgments

This work forms part of the project titled “Estudio autoecológico y modelos de gestión de los rebollares (Quercus pyrenaica Willd.) y normas selvícolas para Pinus pinea L. y Pinus sylvestris L. en Castilla y León”, involving the collaboration of the INIA and the government of Castilla-León.

References

All Tables

All Figures

	Figure 1 Map of sample locations and biogeoclimatic strata (Elena Roselló, 1997).
In the text

	Figure 2 Receiver operating characteristic (ROC) curve for the mortality model (a) and sensitivity-specificity diagram (b).
In the text

	Figure 3 Predicted 10-year probability of individual-tree mortality of pyrenean oak as a function of tree dbh and h/Ho ratio for observed dead and live trees for model data set (a) and (c) and test data set (b) and (d).
In the text

	Figure 4 Relative percentages of observed and predicted 10-year mortality (according to average observed mortality rate, intersection point of sensitivity and specificity, and random number) vs. diameter at breast height (a) and (b) and h/Ho ratio (c) and (d) for model data set (left figures) and test data set (right figures).
In the text