A tree distance-dependent growth and yield model for naturally regenerated pure uneven-aged maritime pine stands in central inland of Portugal

Abstract

• Context

Maritime pine (Pinus pinaster Aiton) is one of the most important Portuguese species, growing in pure stands ranging from even-aged to multi-aged structures. Current growth and yield models were developed only for even-aged, managed stands and/or for very specific regions of Portugal.

• Aims

This paper focuses on the validation of the existing size-class model PBRAVO, adapted to even-aged stands, and on the subsequent development of a single tree distance-dependent growth and yield model (PBIRROL), both in distance-independent and distance-dependent versions, for uneven-aged stands.

• Methods

The new model is composed of four modules, each with a set of sub-models for: tree variable prediction, tree volume prediction, future tree list prediction and growth projection.

• Results

The evaluation of the PBRAVO and PBIRROL models showed that the new model gives more accurate predictions. Moreover, medium-term simulations provided consistent and logical predictions.

• Conclusion

It was verified that individual tree models are more suited to simulate poorly managed uneven-aged stands than diameter distribution models. No clear superiority of distance-dependent models was found over models using just distance-independent measures of inter-tree competition.

1 Introduction

In Portugal, forest land cover represents 39 % of the territory (3.5 × 10⁶ ha) most of which is owned by non-industrial owners (73.4 %). Maritime pine (Pinus pinaster Aiton) is the species that occupies a larger share in the country forest area (28 %). More than half of the Portuguese maritime pine stands (1,003 × 10³ ha) are located in the central region of the country (607 × 10³ ha), where property is very fragmented and where average patch size is less than 5 ha. Portuguese maritime pine usually grows in pure stands (65 %) (681 × 10³ ha), in structures ranging from even-aged to multi aged stands, and 41 % are classified as uneven-aged by the Portuguese National Forest Inventory (412 × 10³ ha). The central region of the country follows the same pattern, with uneven-aged stands corresponding to an area of 259 × 10³ ha. The majority of the cuts applied in these stands are done on the basis of owner’s conjectural economic motivation, rather than on a technical silvicultural prescription (AFN 2010; Alegria 2011b).

Despite pure uneven-aged maritime pine stands having a great expression in Portugal just at the beginning of the 90s of the last century, several studies were carried out in private maritime pine stands of the central inland region of Portugal (Alegria and Tomé 2011). The first studies, conducted in 1992 and 1994, showed that 85 % of these maritime pine stands were unevenly aged—according to the criteria by Páscoa et al. (1981) that classify a stand as uneven-aged if stand tree age variation exceeds 5 years or 10 % of the usual rotation age—due to the fact that they were established by natural regeneration, as well as to the absence of technical management. It was observed that these stands had regenerated over a period of several years, with an average of stand tree age variation of 17 years. In some cases, standard deviations of tree diameters (e.g. dbh, diameter at breast height) were small, and the observed diameter histogram had a typically even-aged shape, despite showing a large within stand tree age variation (Alegria 2011a). This situation had already been reported by several authors who referred that in some cases, uneven-aged stands may have a horizontal structure similar to even-aged stands when they grow on poor sites (Husch et al. 1982; Clutter et al. 1983; Davis and Johnson 1987).

About 71 % of these stands were overstocked, and according to stand stability, an appropriate thinning schedule is needed (Alegria 2011a). Another peculiarity is the way cuttings are applied in these stands. Thinning is applied with no technical rule, often removing the most vigorous tree, which may result in a negative impact for both the potential of residual growing stock and the future stand regeneration. This situation was also reported by Trasobares and Pukkala (2004) in small private woodlots with uneven-aged mixed stands of Pinus sylvestris L. and Pinus nigra Arn. in the northeast of Spain.

These naturally regenerated, over-stocked, stands provide logs with fewer knots, straight and cylindrical bole making them especially suitable for pole production when compared to plantations of this species. This may be considered a good business opportunity. In addition, a study on the economic efficiency of silvicultural scenarios showed that for the existing naturally regenerated maritime pine stands of Portuguese private forest areas, a fully stocked—Wilson’s factor thinning grade D (F _w = 0.23) or E (F _w = 0.25) (Wilson 1946)—pulp wood yield oriented (tree log diameters between 20 and 7 cm) stand prescription should be preferred (Alegria 2011b). In fact, studies on wood properties of maritime pine stands in the central region of Portugal indicate that small-diameter poles have good quality for structural applications, and therefore, future development of national strength grading standards is to be considered. This new utilization may encourage future stand thinning operations by providing extra income to forest owners (Morgado et al. 2009).

To support public and/or private Portuguese maritime pine stand management several stand prescriptions, a set of yield tables and some growth and yield models are available (Alegria 2011b). The existing models were either developed for even-aged, managed stands and/or for very specific regions of Portugal. The first growth and yield model developed for maritime pine in Portugal was a size-class diameter distribution model (PBLEIRIA) for the national maritime pine stands of the Leiria region which are even-aged managed stands. This model was later calibrated for maritime pine stands in Portugal using the data from the 1985/1987 National Forest Inventory and renamed PBRAVO model (Páscoa 1990).This model is the only one that can be applied to maritime pine stands in central inland Portugal. A tree height model and a set of models for tree merchantable volume prediction are also available for the species and region (Alegria and Tomé 2011; Alegria 2011a).

In the view of the importance of naturally regenerated pure uneven-aged stands, it was essential to develop studies on growth and yield dynamics for these specific stands. Therefore, in this study, the PBRAVO model was first validated with data from naturally regenerated uneven-aged stands and checked if reliable predictions could be obtained for this particular management situation (Clutter et al. 1983; Davis and Johnson 1987; Vanclay 1994). Secondly, the existing models for tree height and for tree merchantable volume predictions were calibrated and included as components in a new, more a priori adapted and flexible, single tree distance-dependent growth and yield model, developed both in distance-independent and distance-dependent versions, for the naturally regenerated, pure uneven-aged maritime pine stands of central inland Portugal. Finally, the advantages of using the individual tree model to forecast the development of these stands, instead of the diameter distribution model, were analyzed. A comparison of the distance-independent and distance-dependent options was also explored.

2 Materials and methods

2.1 Data

Naturally regenerated pure uneven-aged maritime pine stands of central inland Portugal were identified based on the data collected from previous studies (Alegria 2011a). Then, stratified sampling, considering strata according to growth variability conditions (Husch et al. 1982; Vanclay 1994), with respect to stand development stage (young—less than 30 years, mature—between 30 to 40 years and old—more than 40 years), site productivity (poor—less than 4.9 m³ ha⁻¹ year⁻¹, medium—between 4.9 to 7.6 m³ ha⁻¹ year⁻¹ and high—more than 7.6 m³ ha⁻¹ year⁻¹) and stand density (low—less than 680 trees ha⁻¹, medium—between 680 to 1,140 trees ha⁻¹ and high—more than 1,140 trees ha⁻¹) was accomplished. Some of the strata considered were not observed in the field, such as very young stands and high productive sites.

From September 1996 to March 1997, 30 semi-permanent circular sample plots of 1,000 m² (radius = 17.84 m) were installed (Fig. 1). Stands unevenness was confirmed if the within stand tree age variability—evaluated by boring a sample of trees—exceeded a 5-year range (10 % of the usual rotation age), as proposed by Páscoa et al. (1981). A within stand tree age variation between 6 and 33 years, averaging 17 years, was observed in the sampled stands. Smaller age variation was observed in young stands (6 to 10 years) while in both mature and old stands, the age variation was greater than 10. Evidence of previous cuts (many stumps in the ground), several small trees (dbh < 5 cm) and some dead trees were observed in almost every sample plot.

Fig. 1

Geographical distribution of semi-permanent sample plots in this study (plots 1999_97) and in previous studies (cutting sites 1993, plots 1994 and plots 1992)

To evaluate growth, recruitment, mortality and harvesting, during the 3 years measurement period, all trees within each sample plot were properly identified by registering their polar coordinates and marked on the field. Sample plots were re-measured from September 1997 to March 1998 and from September 1998 to March 1999.

The following variables were collected on all trees with dbh ≥ 5 cm (2,619) in the sample plots of 1,000 m² (radius = 17.84 m): d, dbh (in centimetres); h, height (in metres); cw, crown width (in metres); hc, height to the crown (in metres) and cr, crown ratio (Table 1). In the first year of measurements (September 1996 to March 1997), a concentric sub-sample plot of 500 m² (radius = 12.62 m) was considered within each plot, and a set of sample trees (314) were selected—two by dbh class and/or a minimum of ten trees per plot—for additional measurements: t, age (years); id5, periodic annual diameter increment in the past 5 years (in centimetres per year) using an increment borer; h _d , bole height (in metres) as merchantable limit and d _h , bole diameter as merchantable limit using the Bitterlich Tele-relascope. The stem diameters (d _h ) at corresponding heights (h _d ) (2,353) were used to evaluate total tree volume (v) and the accumulated merchantable volumes from the ground to each bole section (over bark) (v _m) (Table 1). The total tree volume was obtained by adding stump volume, bole sections volumes and tip volume. Stump volume was evaluated by the cylinder formula; each bole section volume was evaluated with the Smalian’s formula and tip volume with the cone formula (e.g. Loetsch et al. 1973; Husch et al. 1982). From September 1996 to March 1997, the variables d, h, cw, hc, cr, t and id5 were also taken on open-grown trees (52) identified in the study area (Table 2).

Table 1 Summary statistics for the data collected at the tree level in sample plots of 1,000 m² and concentric sub-sample plots of 500 m² in the first year of measurements

Table 2 Summary statistics for the data collected in open-grown trees identified in the study area

Several stand variables were evaluated for the sub-sample plots of 500 m² of area: N, number of trees per hectare (trees per hectare); G, basal area per hectare (square metres per hectare); dg, quadratic mean diameter (in centimetres); $\bar{h}$, mean height (in metres); $\bar{{c}} {r}$, average crown ratio; d _dom, dominant diameter (in centimetres); h _dom, dominant height (in metres); CCF, crown competition factor (in percent); $\bar{t}$, mean age (years) and $\varDelta t$, sample tree age range (e.g. Krajicek et al. 1961; Loetsch et al. 1973; Husch et al. 1982; Clutter et al. 1983; Davis and Johnson 1987; Vanclay 1994) (Table 3).

Table 3 Summary statistics for the stand variables evaluated for the sub-sample plots of 500 m² in the first year of measurements

Study area information regarding elevation, soils, mean annual temperature and mean annual rainfall were also obtained (APA 2007). Elevation ranges between 300 and 1,000 m. Soil types are mainly humic cambisols and eutric lithosols. Based on climatic data collected over a 30-year period (1960–1990), mean annual temperature ranges from 7.5 to 15.0 °C, and mean annual rainfall ranges from 300 to 1,600 mm were observed.

2.2 Evaluation of the PBRAVO model for naturally regenerated uneven-aged stands

The PBRAVO model (Páscoa 1990) uses the Weibull distribution function to predict stand diameter distribution and is composed by a set of equations that make possible to simulate stand projection tables for the following situations: stands previously unthinned, stands previously thinned and stands after the thinning operation. A set of growth equations are also included to project variables in the future. The model includes as input: the average dominant height (h _dom), the age (t) and the number of trees per hectare by diameter at breast height (dbh) class, the diameter classes being organized in a range of 5 cm, with the initial class as [2.5, 7.5[. The output variables of stand table projection are: the number of trees per hectare (N), the basal area per hectare (G), the quadratic mean diameter at breast height (dg), the average height ($\mathop{h}\limits^{\_}$) and the total and merchantable yield by class of timber for industrial use (e.g. round wood, pulp wood and fuel wood). In the simulation process, thinning can be incorporated considering either a systematic thinning or a selective thinning from below (in the dominated stand stratum—trees of smaller dbh, since maritime pine is a light-demanding species). A thinning grade is defined by the stand residual basal area (Páscoa 1990).

The PBRAVO model performance was evaluated using the data collected in this study. Model validation was performed through residuals analysis based on the model error using both graphical visual validation and the calculation of validation statistics. Model bias was evaluated with the mean prediction errors ($\overline{e}$), model precision with the mean of the absolute value of the prediction errors ($a\overline{e}$) and the prediction error variance (${\sigma^2}p$). Modelling efficiency, a measure equivalent to the coefficient of determination, but computed with the prediction errors (R ² P), was also evaluated (Huang et al. 2003). The output variables analysed were: the number of trees per hectare (N) and the diameter distribution, the basal area per hectare (G), the quadratic mean diameter at breast height (dg), the average height ($\mathop{h}\limits^{\_}$), the dominant height (h _dom) and the total volume per hectare (V).

2.3 Model development

The PBIRROL model was developed taking into consideration four modules, each containing a series of sub-models:

• Sub-models for tree variables prediction—tree height, height index for site evaluation, tree age and average tree crown ratio

• Sub-models for tree volume prediction—total volume; tree volume ratio, both to any top height limit and top diameter limit and tree taper compatible to the total volume

• Sub-models for future tree list prediction—recruitment tree list (annual recruitment probability, number of recruitment trees per hectare, recruitment tree diameter and recruitment tree age), mortality tree list (annual tree survival probability) and harvesting tree list (annual tree harvesting probability)

• Sub-models for growth projection—annual tree diameter growth (annual tree diameter potential growth, annual distance-independent tree diameter growth and annual distance-dependent tree diameter growth) and dominant height growth

Modelling was performed using the variables collected and evaluated for the sub-sample plots of 500 m² due to the need to consider a border zone for calculating distance-dependent competition indices. The height index sub-model was developed with the variables collected and evaluated for the sample plots of 1,000 m².

2.3.1 Sub-models for tree variables prediction

The tree height model selection was based on a study for the species in the region (Alegria 2011a). Based on this study, the Prodan (1965) and Harrison et al. (1986) models both modified (Soares and Tomé 2002) were considered as the candidate models for calibration and further analysis in this study.

The classical site indices are based on height–age relationship, but in uneven-aged stands, the height–age relationship can no longer be used to express site productivity because the height growth is not correlated with age and is dependent on the environmental conditions affecting the population during its lifetime (Husch et al. 1982; Huang and Titus 1993; Vanclay 1994). According to McLintock and Bickford (1957), the height–diameter relationship of dominant trees is a reliable approach for uneven-aged stands. In fact, Stout and Shumway (1982) observed that stands of the same species thriving in different conditions have distinct height–diameter relationships. The shape of stand height–diameter curves is steeper in productive sites and flatter in poor sites (Husch et al. 1982). Vanclay and Henry (1988) developed a site productivity measure for uneven-aged coniferous forests in Queensland using a height–diameter relationship that was found to be positively correlated with stand basal area increment with diameter increments of individual trees and, as well as, with several other indicators for site productivity (e.g. stand basal area, periodic annual increment of stand volume and maximum stand height). Field experience suggested that the method worked best in well-stocked monospecific stands (Vanclay 1994). Following this methodology, a height index was tested for site evaluation. The height index was modelled using the data collected in the sample plots of 1,000 m² in the first year of measurements. The guide curve method (Clutter et al. 1983) was used to fit the height–diameter relationship using the monomolecular function modified by Meyer (1940):

$$ h=1.3+A\,\left( {1-{e^{-kd }}} \right) $$

(1)

with A—asymptote and k—species coefficient.

The height index curves were then obtained as

$$ {s_h}25=1.3+\left( {h-1.3} \right)\frac{{1-{e^{-25k }}}}{{1-{e^{-kd }}}} $$

(2)

Sample plots height index for site evaluation S _h 25, defined as the expected tree height at the reference diameter of 25 cm, was computed as the average of the individual tree values s _h 25 (Eq. 2). The use of this height index was evaluated comparing to the mean annual increment of stand volume and the classical site index used in the PBRAVO model.

Models to predict the average crown ratio and tree age were also developed. Tree age was modelled with the tree sample data using a modified version of the logistic function. Later on, average stand age was evaluated and used as an explanatory variable on average crown ratio modelling at the stand level. This was modelled with a modified version of the Weibull function (Soares and Tomé 2001). Further on, this last variable was used as a predictor in the annual diameter increment model.

2.3.2 Sub-models for tree volume prediction

Candidate model selection for tree volume prediction was based on the studies for the species in the region (Alegria and Tomé 2011). The selected models were: the total volume model from Spurr (1952), the volume ratio model to any top height limit from Cao et al. (1980), the volume ratio model to any top diameter limit from Deusen et al. (1981) and the tree taper equation compatible with the total volume equation from Demaerschalk (1973). Thus, these sets of models for volume prediction were calibrated using the variables collected in the set of sample trees (sub-sample plots of 500 m²) (Table 1).

2.3.3 Sub-models for future tree list prediction

To project the actual tree list for the following year, recruitment, mortality and harvesting had to be predicted in order to obtain the future tree list. Recruitment was predicted in a two-stage approach (Vanclay 1994). First, the annual recruitment probability (p(p _I = 1)) was predicted using the logistic function which is limited to the interval of [0, 1] where the dependent variable is a binary variable assuming a value of 0 when the event does not happen and a value of 1 when present.

Then, a conditional linear model was used to predict the number of recruitment trees per hectare (N _I) (Vanclay 1994). Next, recruitment tree diameter at breast height in centimetres (d _I) had to be predicted by Monte Carlo simulation. Finally, recruitment tree age in years (t _I) was fitted using the same model selected for tree age prediction. Mortality was predicted by modelling the annual tree survival probability (p(p _S = 1)) using again the logistic function (Vanclay 1994; Soares and Tomé 2003; Palahí et al. 2003).

In harvesting, silvicultural stand prescriptions for Portuguese maritime pine recommend the application of a thinning from below, since maritime pine is a light-demanding species. But in the stands studied, even though cuts were observed in most of the sample plots, the trees selection for felling did not follow any technical rule. The cuts were correlated with tree dimension especially for the most vigorous trees in the stand. Therefore, a probabilistic approach was used to model current harvesting. Harvesting was modelled in order to predict the annual tree harvesting probability (p(p _C = 1)) using the logistic function (Vanclay 1994). The model can be applied using this harvesting model or another thinning model or rule, in order to simulate alternative thinning types and intensities.

2.3.4 Sub-models for growth projection

The annual diameter increment and the dominant height growth had to be fitted to project growth for the following year. Annual diameter increment was fitted using a potential growth function multiplied by a modifier model expressed as a function of distance-independent and/or distance-dependent competition indices (Vanclay 1994; Peng 2000). The model expressed as a difference equation was the following:

$$ \operatorname{d}{t_2}=\operatorname{d}{t_1}+\mathrm{i}{{\mathrm{d}}_{\mathrm{pot}}} \times \mathrm{modifier}\ \mathrm{function} $$

(3)

First, the model for potential growth (id_pot) was fitted using the data collected on open-grown trees in the study area. The Lundqvist–Korf function (Tomé 2001; Soares and Tomé 2003), expressed as a difference equation with the growth rate parameter as free, was selected to fit potential growth

$$ {i}{{{d}}_{\mathrm{pot}}} = \operatorname{d}{t_2} - \operatorname{d}{t_1} $$

(4)

where $\operatorname{d}{t_2}=A{{\frac{{\operatorname{d}{t_1}}}{A}}^{{{{{\left( {\frac{{{t_1}}}{{{t_2}}}} \right)}}^n}}}}$ with the asymptote $A={b_0}+{b_1}{S_h}25$ and the shape parameter n which influences the moment at which the inflection point occurs.

Subsequently, the annual diameter increment was fitted using a modifier function based on an exponential function of the competition indices (Soares and Tomé 1999) evaluated in this study. Three kinds of modifier functions (MF) were tested: (a) those incorporating distance-independent indices (MFI); (b) those incorporating distance-dependent indices (MFD) and (c) those incorporating both types of competition indices (MFI × MFD).

The distance-independent indices considered were: G, basal area per hectare (in square metres per hectare); G > d, basal area per hectare in trees larger than the subject tree i (in square metres per hectare); cr, tree crown ratio; CCF, crown competition factor (in percent) and several ratios expressing the dimension of the subject tree relative to either the dominant trees or the average trees (e.g. d/d _dom, d/dg, h/h _dom) (Vanclay 1994).

Two types of distance-dependent indices were evaluated: size distance (SD) and point density (PD) (Clutter et al. 1983; Davis and Johnson 1987; Vanclay 1994). Three formulations of the SD type index and one of the PD type index, both in the traditional and the unilateral competition versions, combined with eight search radius criteria for the selection of the competitors among the neighbouring trees j (Tomé and Burkhart 1989; Soares and Tomé 2003), were tested (Table 4).

Table 4 Formulations of the distance-dependent indices, size distance (SD) and point density (PD) in both traditional and unilateral competition versions (_U) and the combined search radius criteria for the selection of the competitors (j) among the neighbouring trees (i)

Finally, dominant height growth was fitted using the Lundqvist–Korf function, expressed as a difference equation with the growth rate parameter as the free parameter

$$ {h_{\mathrm{dom}}}{t_2}=A{{\frac{{{h_{\mathrm{dom}}}{t_1}}}{A}}^{{{{{\left( {\frac{{{t_1}}}{{{t_2}}}} \right)}}^n}}}} $$

(5)

with the asymptote A and the shape parameter n which influence the moment at which the inflection point occurs. This model had been selected for the maritime pine site index curves developed for even-aged stands by Tomé (2001).

2.4 Model fitting and evaluation

Candidate models were evaluated based on best fit, collinearity and prediction performance. Model fitting performance was evaluated with statistics obtained in the SAS output which were: the modelling efficiency (R ²), the adjusted modelling efficiency ($R_{\mathrm{adj}}^2$) and the residual mean square (RMS) (e.g. Myers 1990; SAS Institute Inc. 2004). For models including parameters expressed as linear combinations of the independent variables, the Mallows Cp statistic was also used for model selection (e.g. Myers 1990; Stauffer 2008).

Model multicollinearity analysis was performed based on the computation of variance inflation factors (VIFs) and the condition number of the W′W matrix, where W is the matrix of the partial derivatives in order to each parameter computed at the solution (e.g. Myers 1990; Belsey 1991; Freund and Littell 2000; SAS Institute Inc. 2004). According to Myers (1990), a VIF value greater than 10 was considered to indicate some collinearity in linear models. According to Belsey (1991), a condition number between 5 and 10 indicates that collinearity is not a major problem. In the range of 30 and a maximum of 100, there are problems associated with collinearity, and from 1,000 to 3,000, problems of collinearity are severe. In this study, a condition number greater than 100 was considered as criteria to exclude a model from further analysis.

The analysis of the Press residuals (Myers 1990) was used to evaluate model prediction performance. Based on the Press residuals, the following statistics were computed: the average of the Press residuals (PRESS), the average of the absolute values of the Press residuals (APRESS) and the average of the sum of squares of the Press residuals (SPRESS).

Regression assumptions were checked using the studentized residuals (STR). The assumption of normality of the model error was tested using the Kolmogorov–Smirnov and Anderson–Darling (n > 50) normality tests and by graphical evaluation of the Q–Q plot. Studentized residual plot analysis over the estimated values (e.g. Devore and Peck 1997; Freund and Littell 2000; SAS Institute Inc. 2004) was done to test non-homogeneity for errors distribution variance (Myers 1990). If regression model assumption deviation was found for the best models selected, then the models were fitted again through robust regression techniques, to attempt to overcome the problems of non-normal residuals distribution (Vanclay 1994). To decrease the influence of data points containing large errors in model fit, the iteratively reweighed least squares method (IRLS) with the Huber’s influence function for parameter estimation was used (Myers 1990; SAS Institute Inc. 2004).

In the presence of heterogeneous variance, the weighted regression with a variance estimating function w _i that is inversely related to the variance of errors at each data point was used as well (e.g. Myers 1990; SAS Institute Inc. 2004; Vanclay 1994; Freund and Littell 2000). Binary models, using the logistic function with the stepwise logistic for model selection, were fitted through the maximum likelihood method (e.g. Myers 1990; Vanclay 1994; Freund and Littell 2000; SAS Institute Inc. 2004). Based on the likelihood ratio χ ² test, the most significant explanatory variables were selected. The model performance was assessed using the following statistics: the likelihood ratio test, to analyse the overall significance of the model; the Wald’s test that analyses the significance of the parameter associated to one specific variable when others are present in the model; the odds ratio, calculated for each independent variable and the concordance analysis, based on the analysis of correspondence between real and predicted values which gives an indication of the predictive capacity of the model (e.g. Cody and Smith 1977; Freund and Littell 2000; SAS Institute Inc. 2004).

The definition of the optimal cutoff that will convert the probabilities into (0, 1) binary data was obtained from sorting tables for different levels of cutoff points. Probabilities greater than the cutoff are positive diagnoses to be considered as an event. The statistics of sensitivity and specificity were used for each cutoff: the sensitivity statistic to measure the proportion of true positives that were predicted as events and the specificity statistic to measure the percentage of true negatives that were predicted as non-events. Finally, the optimal cutoff was assessed at the crossing point of sensitivity and specificity curves (e.g. Cody and Smith 1977; SAS Institute Inc. 2004).

The evaluation of the PBIRROL growth and yield model, with the fact that it consists of many sub-models each independently estimated using different techniques, was performed for each individual component and for the system’s final outcome. The tree variables needed as input data for the validation of the PBIRROL model were the dbh (d) and the height of dominant trees, collected in the first year of measurements (September 1996 to March 1997) in the sub-sample plots of 500 m². Sub-model validation was performed through residuals analysis based on the model error using both graphical visual validation and the model bias ($\overline{e}$, $a\overline{e}$, ${\sigma^2}p$ and R ² p) as used for the validation of the PBRAVO model. The output variables analysed were the same as for the PBRAVO model (N and the diameter distribution, G, dg,$\bar{h}$, h _dom and V) and in addition, the average crown ratio ($\bar{{c}} {r}$), the dominant diameter (d _dom) and the mean age ($\bar{t}$).

2.5 Simulations

The PBIRROL model is presented following four modules: sub-models for tree variables prediction, sub-models for tree volume prediction, sub-models for future tree list prediction and sub-models for growth projection. The model structure and operation is synthesised in a flowchart where the sub-model functional relationships are identified (Fig. 2).

Fig. 2

PBIRROL model structure and operation using standard forest inventory data as input to simulate tree and stand variables, future tree list (recruitment, mortality and harvesting) and annual growth (tree diameter and dominant height) to obtain current and future yield on a yearly cycle basis

The model, in its distance-independent version, has the standard data from forest inventory as input: the diameter at breast height of all trees and the height of a subset of sample trees. These data allow stand variable evaluation: number of trees per hectare (N), basal area per hectare (G), quadratic mean diameter at breast height (dg), dominant diameter (d _dom) and dominant height (h _dom). The basal area per hectare in trees larger than the subject tree (G > d) must also be computed for all trees in the sample plot. The model, in its distance-dependent version, also needs the spatial location of all trees in the sample plot as input data (polar or XY-coordinates).

Using the first category of sub-models, the tree height of all trees is simulated allowing the evaluation of the average stand height ($\bar{h}$). The tree variables dbh (d) and height (h) will make site quality evaluation possible. The height index model is applied at the tree level. Then, stand site evaluation is obtained from the average of the individual tree values (S _h 25). The tree age is simulated for stand average age ($\bar{t}$) evaluation, and this last variable is used in average crown ratio simulation. Later on, the stand average age ($\bar{t}$) is also used in the recruitment simulation and the average crown ratio ($\bar{{c}} {r}$) in diameter growth simulation.

The second sub-model category allows stand total and merchantable volume prediction. Tree bole diameters over bark are simulated by using the compatible taper equation to search for the minimum bole height for the dimension class of timber industrial use. Used in conjunction with the volume ratio equation and the merchantable height limits along with the total volume equation, the over bark merchantable volume is simulated. Stand total and merchantable volume are obtained as the sum of individual values and expansion to the hectare.

The third sub-model category is of major importance. Future tree list for the following year will be obtained taking into consideration the incidence of recruitment, mortality and harvest in the stands. First, the probability of having recruitment in a sample plot is first evaluated. If it occurs, then the quantity of trees is assessed. Subsequently, both the dimension and the age of recruitment trees are simulated. These new trees are added to the sample plot tree list. Second, the probability of tree survival is evaluated, to decide which trees will die and have to be excluded from the sample tree list. Third, the probability of a tree being harvested is evaluated. The trees selected for felling will be excluded from the sample tree list.

Finally, the fourth sub-model category is used to project growth for the following year. With the future tree list, annual tree diameter growth is simulated and added to the current tree diameter to obtain the following year tree diameter. Stand dominant height is projected for the following year so that the height of the future tree list can be simulated.

At this point, both the dbh and height of all the trees that will be in the sample plot in the following year have been projected and a new cycle of simulation can be started. Simulations may continue further on, on a yearly cycle basis.

A medium-term simulation, for the possible combinations of stand density (low—N < 680 trees ha⁻¹, medium—N ∈ [680, 1,140[trees ha⁻¹ and high—N ≥ 1,140 trees ha⁻¹) and site quality (low—S _h 25 < 15, medium to high—S _h 25 ≥ 15), as well as with extreme observed conditions, was performed to check if the PBIRROL system’s final outcome predictions were consistency with biological and theoretical behaviour.

Medium-term simulation used the stands selected from the modelling data as a starting point. Using the sample plot data, simulations were started using a 1-year cycle step. A harvesting probability (instead of a thinning rule) as usual in the study area was considered. Simulations stopped when all sample plot trees were finally harvested or dead.

Simulations for low-quality sites used sample plots selected (4) from the modelling data with the following density and age: (a) low stands density—N = 580 trees ha⁻¹ and $\bar{t}$ = 45 years, (b) medium stands density—N = 940 trees ha⁻¹ and $\bar{t}$ = 18 years, (c) high stands density—N = 1,200 trees ha⁻¹ and $\bar{t}$ = 39 years.

Simulations for medium- to high-quality sites used sample plots selected (3) from the modelling data with the following density and age: (a) low stands density—N = 460 trees ha⁻¹ and $\bar{t}$ = 45 years, (b) medium stands density—N = 960 trees ha⁻¹ and $\bar{t}$ = 35 years and (c) high stands density—N = 1,160 trees ha⁻¹ and $\bar{t}$ = 39 years.

3 Results

3.1 Evaluation of the PBRAVO model for naturally regenerated uneven-aged stands

The validation statistics for the PBRAVO model final outcome, on predicting the stand variables N, G, dg, $\bar{h}$, h _dom and V, during the 3-year period when the measurements took place were reasonably good and consistent (Table 5). Modelling efficiency for variables predictions were higher in the first year but started gradually to decrease in the following years, being in the end around 80–87 % for N, dg, h _dom and V. The worst values were found for both G and $\bar{h}$ predictions, 79 and 55 %, accordingly. This situation points to the inability of the PBRAVO model to predict multi-aged stand diameter distribution properly (Table 6), which as a final result affects $\bar{h}$ predictions, and causes the overestimation of stand volume (V). In fact, modelling efficiencies for the number of trees per diameter class were lower than 43 %, showing an overall underestimation tendency for the diameter classes of 7.5 and greater than 27.5 cm, while the classes between 12.5 and 22.5 cm are generally overestimated. This clearly indicates that the Weibull function, used in the PBRAVO model to simulate the diameter distributions, even though it can simulate a normal or even a bi-modal shape diameter distribution, is not able to simulate irregular reverse J-shape diameter distributions as it occurs in the study area.

Table 5 PBRAVO model—validation statistics for the output stand variables simulated using as input data the first year of measurements in sample plots of 1,000 m²

Table 6 PBRAVO model—validation statistics for the diameter distribution simulation using as input data the first year of measurements in sample plots of 1,000 m²

3.2 The PBIRROL model

First, fitting and prediction performances for the set of sub-models that compose each of the four modules of the PBIRROL model are presented. Next, the overall evaluation of the PBIRROL model is analysed. Finally, the medium-term simulations performed with the PBIRROL model are presented.

3.2.1 Sub-models for tree variables prediction

The PBIRROL sub-models selected for tree height, site quality, tree age, average crown ratio prediction (“Appendix 1”) presented good fitting and prediction statistics, but residuals did not follow the normal distribution, with the exception of the average crown ratio model. Therefore, once again, tree height, height guide curve for site quality evaluation and tree age models were fitted using the IRLS method.

3.2.2 Sub-models for tree volume prediction

The selected sub-models for tree volumes prediction (total volume model, volume ratio model to any top height limit, volume ratio model to any top diameter limit and the compatible tree taper model (“Appendix 2”) had both a very good fitting and prediction statistics, but residuals did not follow the normal distribution. The homogeneity of the variance residuals for the case of total volume model was also not observed. Models were fitted again through the IRLS method to overcome this situation.

3.2.3 Sub-models for tree list prediction

Tree list projection models (recruitment, mortality and harvesting; “Appendix 3”) proved to have good fitting and prediction statistics. For instance, the number of concordant pairs for recruitment probability was around 93 %, and for both the annual recruitment per hectare and recruitment tree age models, fitting and prediction statistics were good. The number of concordant pairs obtained for survival probability was around 96 % and for harvesting probability around 76 %.

3.2.4 Sub-models for growth projection

Annual growth models (annual tree diameter increment (over bark) and annual dominant height growth; “Appendix 4”) presented both very good fitting and prediction statistics, but residuals normality was only observed for the annual dominant height growth model. Annual tree diameter increment (over bark), in both distance-independent and distance-dependent versions, was again fitted through the IRLS method. The distance-independent version of the annual tree diameter model uses a modifier function with the following distance-independent competition indices: N, G > d, d/dg and $\bar{{c}} {r}$, while the distance-dependent version of the model uses a second modifier function with the distance-dependent competition index F4H_U. Fitting statistics for the annual tree diameter increment distance-dependent version did not show to be much better than the distance-independent version (“Appendix 4”). Therefore, despite both versions being presented, the distance-independent version was the one considered for further analysis.

3.2.5 Model evaluation

The validation statistics for the PBIRROL sub-models, on predicting the stand variables G, dg, $\bar{h}$, $\bar{{c}} {r}$, d _dom, h _dom and $\bar{t}$, during the 3-year measurements period were good and consistent (Table 7). Good modelling efficiency was observed for N, G, dg and V. Bias for $\bar{h}$, $\bar{{c}} {r}$ and h _dom showed that these variables were slightly overestimated. The performance of the models for $\bar{{c}} {r}$ and $\bar{t}$ predictions was smaller than for the other sub-models. Modelling efficiencies for N, G, dg, $\bar{h}$, d _dom, h _dom and V predictions were higher than 94 % and for $\bar{{c}} {r}$ and $\bar{t}$ around 67–70 %. The diameter distribution showed modelling efficiencies higher than 84 % (Table 8).

Table 7 PBIRROL model—distance-independent version—validation statistics for the stand variables simulated using as input data the first year of measurements in sub-sample plots of 500 m²

Table 8 PBIRROL model—distance-independent version—validation statistics for the diameter distribution simulation using as input data the first year of measurements in sub-sample plots of 500 m²

3.2.6 Simulations

The medium-term simulations performed showed that the PBIRROL model has a consistent biological behaviour. To exemplify the model’s behaviour, simulations for stand density (N), total volume (V) and quadratic mean diameter (dg) are presented (combinations for low-quality site and low, medium and high stand density) (Fig. 3). The simulations ended when all sample plot trees were harvested or dead, which happened at the average age of 60, 39 and 49 according to low, medium and high stand densities. Stand densities have decreased in time, despite now and then showing some small increase due to recruitment. Stand volume has increased in time, now and then showing some small decrease due to mortality and/or harvesting. Stand average diameter has also increased in time, showing now and then some small decrease due to the selection of larger and/or older trees for harvesting, which is consistent with what is observed in this type of stands.

Fig. 3

Medium-term simulation for the variables stand density (N), total volume (V) and quadratic mean diameter (dg) for the stands selected from the modelling data with low-quality site (S _h 25 < 15 m) and low stand density (N = 580 trees ha⁻¹), medium stand density (N = 940 trees ha⁻¹) and high stand density (N = 1,200 trees ha⁻¹)

4 Discussion

The validation of the existing size-class PBRAVO model, in comparison with the evaluation of the PBIRROL model, showed that the PBIRROL model allows more accurate predictions for the naturally regenerated, pure uneven-aged maritime pine stands, in central inland Portugal. In fact, higher modelling efficiencies for the stand variables N, G, dg, $\bar{h}$, h _dom and V were observed using the PBIRROL model. Moreover, it was clearly observed that the major weaknesses found on the PBRAVO model were both the inability to properly predict multi-aged stand diameter distributions and to overestimate stand volume. Although the Weibull distribution function used to predict stand diameter distribution has the ability to simulate normal, bi-modal and irregular shape diameter distribution, the PBRAVO model could not successfully predict the reverse J-shape diameter distributions as it occurs in the study area. This situation has a direct impact on stand variables prediction, such as $\bar{h}$ and V, since they are evaluated using the simulated diameter distribution data.

It should be noted that in some size-class (diameter) distribution models developed for uneven-aged stands, the reverse J-shape distribution (well suited to sustainably manage shade-tolerant species) (Davis and Johnson 1987) has been used; another approach in use was the stand projection table by Ek (1974) that consists of equations to predict periodic recruitment, mortality and survivor growth by 5-cm-diameter classes (Peng 2000).

Even though size-class models have the advantage of requiring only overall stand values as input and provide some detailed size class information, they are not flexible enough to evaluate a broad range of stand treatments. By contrast, individual tree models provide maximum detail and flexibility that allows alternative utilization and stand treatment evaluation. For structurally complex forests, such as multi-aged structures and/or mixed stands, individual tree models are more flexible and reliable than stand models or size-class models. Moreover, a number of growth and yield models, using neither age as variable or site index for assessing site quality, have been proposed for uneven-aged stands in USA and in Austria. These modelling techniques can also be applied to even-aged stands (Peng 2000).

The PBIRROL model combines some of the approaches in use in growth and yield models for uneven-aged stands. As an example, the PROGNOSIS model (Stage 1973; Wykoff et al. 1982), one of the few individual tree distance-independent models that does not need site index and age as variables, is able to predict growth and yield of forest stands of any composition, from pure even-aged to mixed uneven-aged structures, in North Idaho, USA. Its growth equations are a function of tree size, vigour and dominance; site productivity is described by site characteristics (Peng 2000). Another example is the FOREST model (Ek and Monserud 1974), an individual tree distance-dependent model designed to simulate the growth and reproduction of even or uneven-aged mixed species stands of northern hardwoods in USA. In this case, spatial information is summarized by a competition index, which is calculated for every tree at the start of each growth period to account for availability of light, moisture and nutrients to the tree (Peng 2000).

The PBIRROL model uses a height index model for site quality evaluation that worked very well, being positively correlated with both mean annual increment of stand volume (r = 0.71) and the classical site index used in the PBRAVO model (SI₅₀) (r = 0.84). The height index, developed in this study, can be applied without the need for age data, which is considered to be an advantage when compared with the use of a classical site index.

For tree volume prediction, a tree volume model, a tree volume ratio to any top height limit model, a tree volume ratio to any top diameter limit model and a compatible tree taper model were proposed. These models allow total and merchantable volumes predictions, for different definitions of wood for industrial use (e.g. round wood, pulp wood and fuel wood) (Alegria and Tomé 2011). It is worthwhile to note that the PBRAVO model uses the data organized in a diameter frequency table to predict total volume with a tree volume model and the merchantable volumes, for different definitions of wood for industrial, by integrating a tree taper model with diameter class-specific parameters (Páscoa 1990).

Future tree list prediction considered the components of forest dynamics: recruitment, natural mortality and harvesting. The approach used for recruitment and mortality is widely used (Trasobares et al. 2004; Pukkala et al. 2009), but the same cannot be said for the harvesting component. A probabilistic approach was used to solve the fact that these stands are conducted based on the owner’s conjectural economic motivation, rather than on a technical silvicultural prescription. However, the model can be used with different thinning types and intensities which will allow for an economic analysis of the advantages/disadvantages of the present thinning method, as compared with alternative ones, namely the low thinning which is common in national even-aged Portuguese maritime pine stands management.

Finally, for growth projection, two growth models were considered, one for tree diameter growth and another for the growth of the dominant trees. Annual tree diameter growth, in both distance-independent and distance-dependent versions, was modelled. Fitting and predictions statistics for both versions pointed out to a small improvement for the distance-dependent version. These results are consistent with other studies (Vanclay 1994; Burkhart and Tomé 2012) suggesting that the distance-dependent competition indices rarely provide better estimates of diameter increment than some simple measurements, such as the basal area per hectare of trees larger than the subject tree i (G > d). In fact, this last variable was included in several sub-models of the PBIRROL model. Furthermore, the model in its spatial version would be very limited in use, since no tree coordinates are collected in standard forest inventory in Portugal, unless a tree stand structure simulator was available.

However, some limitations on data collected should also be addressed: the short period of passive stands monitoring and the lack of data from experimental plots, which were inexistent. Moreover, the measurement period of 1 year might not have always been enough to guarantee that the growth exceeded the measurement errors, as recommended by Vanclay (1994).

Despite of this, the PBIRROL sub-models exhibited good fitting and prediction statistics and a logic biological behaviour. The validation results for the stand variables G, dg, $\bar{h}$, $\bar{{c}} {r}$, d _dom, h _dom and $\bar{t}$, during the 3-year measurements period, also proved to be good. Medium-term model simulations showed that consistent and logical predictions are provided according to existing forest growth and yield theory (e.g. Loetsch et al. 1973; Clutter et al. 1983; Davis and Johnson 1987).

The PBIRROL model validation, using an independent validation data set, should be undertaken as growth data are collected and sub-models improved under changing conditions or operating environment. Third party validation should be also promoted (Huang et al. 2003). However, this has not yet been possible due to the inexistence of more growth data for this type of stands.

Finally, the PBIRROL model and/or some of its sub-models are expected to have an important role to support the implementation of the forest management plans required by law. Besides, the methodology developed here for maritime pine in Portugal can be used for other species and regions in which naturally regenerated uneven-aged stands are important.

Appendices

Appendix 1

Table 9 PBIRROL model—sub-models for tree variables prediction (tree height, height index, tree age and average tree crown ratio)

Appendix 2

Table 10 PBIRROL model—sub-models for tree volume prediction over bark (tree volume, tree volume ratio to any top merchantable limit and tree bole diameters)

Appendix 3

Table 11 PBIRROL model—sub-models for future tree list prediction (recruitment, mortality and harvesting)

Appendix 4

Table 12 PBIRROL model—sub-models for growth projection (annual tree diameter growth in both distance-dependent and distance-independent versions and dominant height growth)