Abstract
While effects of thinning and natural disturbances on stand density play a central role for forest growth, their representation in large-scale studies is restricted by both model and data availability. Here a forest growth model was combined with a newly developed generic thinning model to estimate stand density and site productivity based on widely available inventory data (tree species, age class, volume, and increment). The combined model successfully coupled biomass, increment, and stand closure (=stand density/self-thinning limited stand density), as indicated by cross-validation against European-wide inventory data. The improvement in model performance attained by including variable stand closure among age cohorts compared to a fixed closure suggests that stand closure is an important parameter for accurate forest growth modeling also at large scales.
1. Introduction
A common challenge in large-scale forest planning is to find an optimal balance among different land use management options, such as forest carbon sequestration, timber production, and other forest-based services. In such an analysis, regional-to-national and even continental-scale computationally efficient modeling of potential effects of forestry activities and disturbances on these services is necessary. These effects depend on many factors of which forest productivity, stand density, and characteristics of management, such as thinning and rotation length, are particularly important [1, 2]. In addition, stand density and thinning are important in relation to other forestry topics, such as carbon sequestration in forests [3], biodiversity [4], sensitivity to wind damage [5], sensitivity to insect damage [6], and fire risk [7].
Stand density and thinning have received considerable attention in forest research, and a large body of knowledge has been accumulated on the effects of thinning on physiology and productivity from the leaf to the stand level [8]. However, despite the central role of thinning in forestry, large-scale (continental) forest scenario analysis tools hardly incorporate state-of-the-art knowledge of thinning and stand density effects [9–11]. A reason for not including these effects in large-scale analyses may be the site and species-specific nature of most stand-scale stand density models. The stand scale models are often limited in their generality and require input data that limits their applicability over large scales where coherent data are scarce. Stand and individual-based thinning models relying on empirical findings of growth and yield studies, such as the MELA system [12], Motti [13], and Prognaus [14], have to date not been applied and validated at continental scale, to our knowledge.
Even when appropriate models of thinning and stand density effects are available, the models generally require data on stand density or thinning intensity for initialization. These data are often not available at continental scales. For example, European-scale forest inventories [15] have been widely used for forest simulations, particularly using the EFISCEN model, for example, [16, 17]. These inventories do not contain information on stand density or thinning regime, which in the case of EFISCEN were derived from local handbooks and expert knowledge [18]. Such methods of deriving stand density and thinning regime have drawbacks in terms of accessibility and consistency of input information over large areas. Stand density over large areas can also be estimated from plot data of individual tree sizes, requiring detailed data sets [19] not available over large scales. Remotely sensed data types do not always yield reliable results (e.g., applying conventional optical remote sensing methods [20]) or are technically and financially demanding, for example, using airborne lidar [21, 22]. In summary, the availability of stand density data over large scales is limited by a number of factors. Thus, an indirect method of estimating stand density based on readily available data could significantly improve forest characterization and subsequent simulation of growth and management.
Here, a method is suggested for obtaining site productivity and relative stand density (closure) from inventory data on stand volume and increment. The underlying hypothesis is that there is a causal relationship among age, stand volume, increment, stand density, and site productivity. This relationship is described by a model that combines a growth model for unthinned stands with a simple mechanistic model of the thinning effects (closure effects) on growth and mortality, recently developed [8]. Due to the use of a relative measure of stand closure (closure = density/self-thinning limited density) in the thinning effects model, effects of site productivity and age are factored out, making the thinning effects model applicable to all ages and sites. The effects of site productivity and species-specific growth patterns are added through the growth submodel for unthinned stands, which is parameterized for different species using yield table data. Based on the combined model, relationships between biomass and increment in inventory data can be used to infer stand density from these widely available data by inverse modeling. The goal is to obtain a large-scale modeling framework for prediction of stand density and site productivity based on inventory data. Estimated data on stand density and productivity can subsequently be used to initialize simulations of forest growth and management scenarios, including effects of thinning and stand density over large scales. European-scale forest inventories [15] were used to test the model.
2. Materials and Methods
2.1. Model Structure Overview
The large-scale model consists of two mechanistic-deterministic submodels: the growth model and the thinning effects model. These submodels are parameterized and evaluated separately based on separate datasets before they are combined in the final large-scale model. The growth submodel predicts growth of unthinned stands in response to a site productivity variable and a number of fixed species-specific parameters, which were estimated based on yield table data for each species. The thinning submodel [8] defines how growth and mortality are affected by stand closure (), which is equal to the current stand density divided by the maximum possible stand density (the self-thinning limit) at a fixed tree size. Because thinning is modeled in terms of relative effects on growth and mortality it can be combined with the growth submodel by multiplication. The combined model predicts biomass and increment in response to the variables site productivity () and closure (), while other parameters (including the species-specific parameters) remain fixed. This model was evaluated at large scale by testing its ability to predict observed increments and volumes as reported in a European-wide forest inventory dataset. Variables and parameters are described in Table 1.
2.2. Growth Submodel for Unthinned Stands
The growth submodel is based on equations for stem wood growth because this is generally the most interesting part of the tree for forestry and because it is a good predictor of total tree biomass. The starting point is a differential equation for growth of stem biomass per hectare () similar to the Richards equation [24] with modifications to allow the growth peak to shift on the biomass axis as a function of site productivity and to link productivity to asymptotic (maximum) biomass () where and is the productivity index (unitless) of the site, which integrates effects of climate and soil productivity on growth and which affects both the initial increment (1) and (2) of the stand. is a scaling factor, which enables to be changed without altering the maximum growth rate. controls the shape of the growth rate-biomass curve, that is, the biomass at maximum growth rate relative to the maximum biomass. The constant exponent controls the effect of on growth rate. For the species studied here, estimated was positive, which means that the maximal growth rate increases faster than linearly with and that peak growth rate occurs earlier the more productive a site is. and are species-dependent constants that control and the initial growth rate relative to p, respectively. An example of the effects of differing values of and between two species with contrasting growth curves is illustrated in Figure 1.
Equation (2) means that is proportional to site productivity () for a given species in the absence of thinning (or other disturbance). In the absence of data suitable for empirical testing of this relationship, its plausibility was tested using a physiological model [25], which produced approximately linear relationships between maximum sapwood mass and wood increment in response to soil N availability and light availability (Figure 2). Thus, (2) is a simplified yet plausible first-order approximation.
Solving (1) for gives the biomass as a function of age (). In (3) is the initial at , and the subscript indicates a closed (unthinned) stand.
This stand growth model is only valid for even-aged, unthinned (closed) forest stands. Thus, before the stand growth model can be analyzed and parameterized for managed (thinned) stands, the effects of mortality and stand density must be added to the framework.
2.3. Self-Thinning and Mortality
Because self-thinning, mortality, and thinning effects are presented in detail in [8] only a brief description is given here. To model the forest stand density during the development of an unthinned stand, a self-thinning equation predicts maximum density(4) In (4) is the number of trees, and is the biomass of the average tree. For a fixed area, the number of trees () is decreasing as and the total stand biomass () are increasing. Although the self-thinning exponent (, (4)) is not universal among species and site conditions, there were no statistically robust and site-independent species differences for the majority of species in this study [26]. Thus, a single value of was used for all species [8]. The self-thinning parameter determines the absolute value of for a given , however, because in this study only relative density is used and never absolute numbers of have no influence on our results.
As shown in [8] mortality due to self-thinning in a closed stand (, (5)) is proportional to the increment of live biomass () In (5), is the size of self-thinning trees relative to the size of the average tree. It is natural that is less than 1, since smaller trees are suppressed by larger competitors due to size-asymmetric competition for light and therefore are more likely to die in the self-thinning process than larger trees.
2.4. Thinning Effects
On the stand scale, a thinning response can be divided into two effects. First, the total stand biomass and production (NPP) is reduced due to reduced total growing stock and resource capture [27]. Second, the mortality due to competition (self-thinning) is reduced, which is strongly linked to the improved growth of the remaining trees [28]. To model these effects in a way that will be valid across different sites and productivities, the concept of relative density or closure (; [29]) was used, that is, the number of trees relative to the maximum number of trees (given by the self-thinning limit) for a fixed individual average tree size (). If the removed trees are of a different mean size than the remaining trees, of the thinned stand can be calculated based on the new mean tree size after the thinning (see [8]):
The reduction in stand growth rate through reduced resource absorption is related to the reduction in closure () according to (6), where is the ratio growth at closure to growth at [8]. In addition to the growth rate, also the self-thinning rate is affected by . As further discussed in [8], self-thinning mortality declines steeply with declining due to decreasing competition among trees. This effect is modeled in terms of self-thinning of an unthinned stand relative to self-thinning of a closed stand () given by (7), which captured well empirical observations of different stands and species subject to thinning experiments [8] Mortality other than self-thinning is included in the growth model as density-independent mortality, see below.
2.5. Growth and Mortality Dynamics of Thinned Stands
The total effect of a thinning event on net increment () depends on the relative strength of the effects on stand growth and mortality [8]. In an unthinned (closed) stand, net increment () and self-thinning mortality () are linked according to (5). In a thinned (open) stand, net increment is a function of gross production (8), self-thinning (9), and density-independent mortality (; (10)). Compared to a closed stand, gross production is reduced with according to (6), and self thinning is reduced according to (7) in a thinned stand Density-independent mortality (, (10)) has been shown to be higher for very small and very big trees [30], that is, during slow growing stages (in terms of net growth rate) of stand development and is therefore modeled to be a function of the net increment () relative to its maximum over the stand development () A constant maximum density-independent mortality % per year across species and countries was assumed. In general, this is an unrealistically simplistic assumption. However, given the focus of this study on managed forests where density-independent mortality is small due to the absence of very old stands and management practices to avoid mortality in young stands, the value of (within realistic bounds) has little impact on the results. In (10), is obtained by maximizing (1) with respect to . By combining (9) and (10), equations for net (11) and gross increment (12) for thinned stands were derived. In these equations is given by (1). To interpret time-series type data, such as the yield tables used here, it is not sufficient to know the effects of a single value on growth and mortality but it is also necessary to predict the development of over time. In [8], a thinning scenario was presented that aimed at maximizing net biomass production by keeping at an optimal average value by successive thinnings. In that scenario, is reduced by thinning to and thereafter increases again to eventually reach a value , where a new thinning is triggered. The mean over time during such a scenario with one or more thinning will be approximately equal to .
The stand biomass for a single thinned or disturbed stand as a function of time, , will be highly variable depending on whether it is observed shortly after a thinning or if the stand is nearly fully closed (). However, for averaged over many stands (such as for the inventory data used here) subject to reasonable thinning scenarios (<50% volume removal, ) in managed forest (not old growth), can be approximated by of a unthinned forest (3) times the temporal mean of (13) if the thinning scenario is unknown In the inventory data that we use to test the model, all stands of the same species, of the same age, and in the same region are aggregated and modeled by a single and by either a single or age-dependent s. While in reality there is heterogeneity among such stands, which is not captured by our approach, our hypothesis is that the mean variables used in our model: closure () and site productivity (representing climate and soil factors) are reasonable predictors of mean biomass and increments. The linkages between the model’s variables, which are subjected to the testing based on the inventory data, are illustrated in Figure 3.
2.6. Model Parameterization
2.6.1. Thinning Parameters
Using the hybrid patch model PICUS v1.41 [31] which simulates mortality at the individual tree level, the value of the parameter describing the size of self-thinned trees relative to surviving trees in an even-aged stand () was estimated to vary between 0.5 and 0.25 in a generic simulation experiment [32]. Because an empirical study showed similar or slightly lower values [8] and because there is not sufficient data available for species-specific estimation, was chosen in this study.
The effects of thinning on growth and mortality are determined by the parameters and , which have been parameterized in a previous study (Table 1, [8]), using experimental data from the literature. In the absence of data for species-specific parameterization of and for the species modeled here, common values of and for all species were determined from experimental data [8]. Fortunately, experimental results [8] suggest that the effects of variability among species in and on growth and mortality are relatively small, particularly in comparison to the effects of other species differences, for example, in growth rate and maximum attainable biomass.
2.6.2. Growth Parameters
The parameters of the growth functions are estimated for ten tree species by fitting yield table data to gross and net increment (12) and (11). The yield table data summarized in [23] is based on empirical information on age, stem volume, stem volume increment, and stand density () at different site productivities for each species, and were published between 1950–2000. For each species, yield tables representing a range of site productivities for all European countries available (21 countries) were combined to make the estimated species parameters as representative as possible. Because the model was parameterized independently for each species using species-specific yield table data, the value estimated for one species at a given site may not be valid for another species at the same site. However, this restriction has no effect on the results presented here because no species substitutions were simulated.
The yield table volumes and increments were converted to biomass () using wood density (Table 2). As the yield tables represent lightly (<10% of volume) but frequently thinned stands, where potential self-thinning most likely was avoided, a maximum closure, , was used rather than the used in the optimal scenario for stronger nonselective thinnings obtained in [8]. A higher leads to higher biomass production if self-thinning mortality can be avoided. The yield table data on net increment () and gross increment () were fitted to (11) and (12) to estimate the species parameters () for each species while other parameters were kept fixed. A common was used for all species. The parameters of the growth functions (Table 2) were estimated using a Bayesian model with standard noninformative priors on parameters [33], and the program was run on WinBUGS [34]. Measurement error was assumed to be a zero mean Gaussian. While the site productivity index () is allowed to vary between stands (yield tables), the species constant parameters are assumed to be independent of any other factors, such as region, site, or climate. Initial biomass () was calculated from species-specific stand density recommendations for Central Europe [35]. An example of the parameter fitting and the effects of differences between species is given in Figure 1.
2.7. Large-Scale Model Testing
2.7.1. Dataset
The combined model (combining the species-specific growth functions and the thinning effects model) was tested using a European-wide forest inventory (EFI database [15]), which contains standing volume, increment and area for each forest cohort (defined by age class (e.g., 0–10 yrs, 10–20 yrs, etc.), species (), and region (, 1–22 regions per country, 20 countries; Table 3). For the analysis the data was organized as age-structured forest units, which is defined as the group of all age cohorts belonging to the same species and region. Thus, for each forest unit the data consists of biomass and increment for each age class, used here as a time (age) series of biomass and increment data (a space for time substitution). A summary of the data is given in Table 3.
2.7.2. Model Parameterization
The model was fitted to the inventory data by adjusting closure () and site productivity index (), while all species-specific growth parameters and the parameters for thinning effects were fixed. The productivity index was assumed not to vary among age cohorts for the same species and region (within a forest unit). For , two alternative model variants were compared to test if variable among age cohorts (model B) improved the model compared to a fixed (model A). The estimation of and for each forest unit was done by fitting data on gross increment () versus biomass () to (12), and versus age () to (13) by minimizing mean squared error of modeled versus measured values (using MathCad software). In the error function, increment was weighted by a factor 0.25 and biomass by 1 to account for the larger relative uncertainty in increment than in biomass (volume) measurements. The best fit parameters were calculated using a grid search algorithm that always finds the global minimum of the error function (with a 1.5% error margin in and ), avoiding potential problems of local minima.
To avoid strong effects of differences in initial density (), for which measured data is lacking, cohorts younger than a minimum-age were excluded from the analysis. This minimum age was species specific and was calculated as the age when gross increment () of an unthinned stand reaches half its maximum value (based on (1) and (3)). A lower closure limit of (very low but possible for young sparsely planted stands) was set to avoid prediction of unrealistically low closure.
2.7.3. Model Evaluation
Cross-validation was used to evaluate the performance of the models, which is a widely used method for estimating prediction error. It allows comparison of completely different models and is independent of the number of parameters and possible correlation between them as well as of the distributional assumptions [36]. For each forest unit, data points for biomass and increment () were dropped one at a time while the remaining points (the other age classes) were used to estimate the and parameters (see Section 2.7.2). The estimated and values were then used to predict biomass and increment of the dropped data points (that were not used in the parameter estimation). The differences between these predictions (of the dropped data points) and the data were used to calculate percentage root mean square error (PRMSE), which was used to evaluate the model’s ability to predict biomass and increment, averaged for all data (Figure 5) for each species (Table 4), and for each country (Table 3).
3. Results
3.1. Growth Submodel Parameterization
Modeled net and gross increment (11) and (12) fitted to data from yield tables was used to estimate species-specific parameters of the growth submodel for unthinned stands (Table 2). The model fit to data was generally better for net increment ( to 0.93) than gross increment ( to 0.82). There was a tendency of the model to underestimate at high gross productivity, which may be due to differences in the thinning regime in terms of the maximum closure of the stand (), which, in the absence of data, was assumed to be constant. The results indicate that the model based on species-specific parameters combined with a single site productivity parameter describes well net increment of the tested species as given by yield tables, whereas the prediction of gross increment at high biomass is reasonable but less accurate than for net increment (Figure 4).
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
3.2. Large-Scale Model Testing
The model was fitted to measured biomass, increment, and age data for a total of 4417 forest cohorts in 432 forest units (delineated by region and species), first with a common closure for all age cohorts in a forest unit (fixed , model A) and second with individual values of for each age cohort (model B). The results in Figure 5 show a large improvement in model-data agreement using variable compared to fixed . For biomass estimates, cross-validated percent root mean squared error, PRMSE = 27 and 22, for model A and B, respectively, and for increment PRMSE = 38 and 30, for model A and B, respectively. Model B also shows lower mean bias in fitted biomass and increment than model A (Figure 6). Thus, as indicated by the cross-validation, allowing to vary with age improves the model in terms of its representation of the relationships between age, biomass, increment, and closure compared to using a fix value of . Due to the better performance of model B than model A the following results and discussion will be focused mainly on model B.
(a)
(b)
(c)
(d)
The model including age data and variable (model B) performed better than model A for all species, with cross validated PRMSE for biomass between 15.5 for Abies sp. and 27.5 for Pseudotsuga (Table 4). Modelled (model B) mean closure across all regions and ages varied among species between 0.52 for Populus sp. to 0.76 for Betula pendula.
Model performance for model B was similar among countries (Table 3), with the best results for Denmark (PRMSE = 9.8 and 11.9 for biomass and increment, resp.) and the worst for Finland (PRMSE = 30.5 and 55.8 for biomass and increment, resp.). Comparing mean closure () among countries reveals a tendency for higher values in eastern Europe than in western Europe (Table 3).
4. Discussion
4.1. Growth Submodel
The growth submodel for closed stands was able to reasonably explain growth for most species in terms of net and gross increment although values for gross increments in Abies sp. and Picea abies were low (Table 2). Because the yield table data represent empirical data from a wide range of site productivities for each species (e.g., due to variation in soil fertility and climate), the estimated fixed growth parameters should be representative for most European forests. Most likely, this growth modeling could be further improved by allowing other parameters in addition to to vary in response to site conditions. However, without additional site data this would only serve to increase uncertainty in the interpretation of stand density estimates as discussed below. Because focus of this study was the use of the growth submodel in combination with the thinning model (to estimate stand density) rather than as stand-alone application, its merits are discussed further in terms of its effects on the main results, that is, the productivity and closure estimates (see below). The thinning submodel has been thoroughly discussed in an earlier paper [8].
4.2. Large Scale Model Performance
The improvement in the agreement (cross-validated PRSME) between measured and modeled data obtained by including age-related variation in closure (Figure 5) indicates that increment and biomass are strongly dependent on stand closure (), as represented in the model. This conclusion is further supported by the result that mean estimated follows the expected development over time, increasing as stands close after planting and decreasing in response to onset of thinnings and increased risk for disturbances at higher age (Figure 7) [37]. It could be argued that the improvement in model performance obtained by including variable could instead have been achieved by letting productivity () vary among age cohorts. This is however not possible because in the model, affects increment and biomass almost equally (see Section 2) and can therefore not explain the variation in the measured biomass versus increment that caused the lower PRMSE of model A (Figures 5(a) and 5(b)). In addition, there is no reason to expect to vary significantly among age cohorts, each representing means over many forest sites, except for old forests as discussed below and in case site productivity has changed systematically over time. Such systematic trends, for example, related to nitrogen deposition, atmospheric CO2, and temperature, are outside the scope of the current analysis but are interesting topics for further research. Although many temperate/boreal species may respond similarly to thinning [8], it is likely that the model could be improved by including species-specific thinning/stand density effects on growth and mortality. However the data necessary for such a species-specific parameterization of thinning effects is currently not available [8].
4.3. Age-Related Biases
Despite the overall good agreement between modeled and measured for model B there are slight systematic deviations. There is a tendency of the model to underestimate at low values and young age and overestimate at high values (Figure 6(d)). At the same time, predicted increment shows a small but opposite pattern of bias, with slight over estimation at young ages and underestimation at high age. This pattern of biases suggests that for some young or old forests the model cannot reconcile the observed relationship between biomass and increment. A likely factor in the underestimation at young age is a too low initial biomass (). For example, there may be more trees per hectare in naturally regenerated stands compared to planted stands, where the latter was assumed in the current study. In addition, as it was observed in many cases that the measured yearly increments for young ages were not sufficiently high to explain (by multiplication with age) the inventory biomass, measurement irregularities affecting the inventory data may have contributed to these disagreements between model and data for young forests.
A potential source of bias is that predicted biomass by model B may be affected by systematic difference in site productivity () among age classes. As forests at sites with very low productivity grow slowly and are harvested later than at high productivity sites, a cohort that is older than the median rotation age is likely to come from a site with lower than the average for the region. Thus, variability in will lead to a distribution over age that is declining at ages older than the mean rotation age. This lower at old age leads to overestimation of modeled at old ages and underestimation at young ages, since a constant was assumed for all ages.
4.4. Differences among Countries
Although the variation in closure among countries should be interpreted with care due to the variable number of forest units among countries, the results indicate that mean closure is higher in eastern than in western European countries (Table 3). The high mean estimated for Portugal is an exception to this trend, which however is based on a single forest unit with three age classes. The trend may be related to the large cohorts of regrowing forests, following large-scale afforestations in eastern and central Europe after World War II, that have not been thinned due to a lack of industrial capacity to process the thinned trees as well as socioeconomic aspects [38]. The lowest mean closure () was estimated for Finland which may indicate a very intensive thinning regime or disturbance effects not included in the model (see above, Section 2.5).
4.5. Limitations and Applications of the Model
The approach necessarily sacrifices complexity and details compared to comprehensive individual-based models, for example, [14, 39, 40] and the limited species (e.g., non-species-specific thinning responses), and site information used in the model makes the model unsuitable for practical, stand-level forestry applications. However, the selected mean tree approach is a stand-structural advancement of the state-of-the-art compared to structurally simple scenario tools applied at continental scale (e.g., [18]). In addition, the model is unique in its representation of thinning and stand density effects and in the range of countries and species for which it has been tested. Due to its simplicity and generality, the thinning effects submodel can be used to add effects of stand density and thinning to a wide range of existing forest growth models, as discussed in [8]. As shown here, the resulting combination of thinning and growth models can be used to obtain initialization data on stand density based on widely available inventory data on volume (or biomass) and increment. In addition, such models could be used with remotely sensed data types for productivity, for example, [41] and biomass, for example, [42]. The estimated stand density data can then be used to improve predictions of future forest development in comparison to models that do not account for the effects of stand density. Furthermore, a model including stand density effect can be used to simulate the effect of different thinning scenarios in order to evaluate the management options in large-scale land-use models.
The model has been used to calculate potentials for timber production and carbon storage in a European-wide land use optimization study [43]. In this work the model (model B) was used to initialize biomass, increment, and relative stand density for all forest units. Thereafter the development for each forest unit under different thinning scenarios was simulated.
5. Conclusions
A site-independent model of closure effects was combined with a growth model for unthinned stands that accounts for the effect of species differences and site productivity. Evaluation against European-wide inventory data showed that the resulting combined model can reasonably predict the relationships among biomass, increment, age, and closure over large scales in managed forests. The large improvement in model agreement with data attained by including variable closure among age cohorts compared to a fixed closure suggests that closure is an important parameter for accurate forest growth modeling, also at continental scale. Consequently, large-scale forestry models could be significantly improved by including stand density (closure) effects. Due to the simplicity and the generality of the formulations of stand density effects (i.e., site and age independent), the approach presented here can readily be used to add thinning and stand density effects to other forest models.
Acknowledgments
The authors thank Guenther Seufert at JRC and Gert-Jan Nabuurs at Alterra and EFI for their assistance in making the yield table and forest inventory datasets accessible to us. This research received funding from the European Community’s 7th Framework Programme (FP7) under the Grant no. 212535, Climate Change-Terrestrial Adaptation and Mitigation in Europe (CC-TAME), http://www.cctame.eu/, and from the 6th Framework Programme (FP6) under the Grant no. SSPI551 CT-2003/503614 (INSEA).