Growth and yield model for uneven-aged mixtures of Pinus sylvestris L. and Pinus nigra Arn. in Catalonia, north-east Spain

Annals of Forest Science, Jul 2018

Antoni Trasobares, Timo Pukkala, Jari Miina

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

http://www.afs-journal.org/articles/forest/pdf/2004/01/F4102.pdf

Growth and yield model for uneven-aged mixtures of Pinus sylvestris L. and Pinus nigra Arn. in Catalonia, north-east Spain

Ann. For. Sci. Growth and yield model for uneven-aged mixtures of Pinus sylvestris L. and Pinus nigra Arn. in Catalonia, north-east Spain Antoni TRASOBARES 2 Timo PUKKALA 0 ri MIINA 1 0 University of Joensuu, Faculty of Forestry , PO Box 111, 80101 Joensuu , Finland 1 Finnish Forest Research Institute, Joensuu Research Centre , PO Box 68, 80101 Joensuu , Finland 2 Centre Tecnològic Forestal de Catalunya , Pujada del seminari s/n, 25280, Solsona , Spain - A distance-independent diameter growth model, a static height model, an ingrowth model and a survival model for uneven-aged mixtures of Pinus sylvestris L. and Pinus nigra Arn. in Catalonia (north-east Spain) were developed. Separate models were developed for P. sylvestris and P. nigra. These models enable stand development to be simulated on an individual tree basis. The models are based on 922 permanent sample plots established in 1989 and 1990 and remeasured in 2000 and 2001 by the Spanish National Forest Inventory. The diameter growth models are based on 8058 and 5695 observations, the height models on 8173 and 5721 observations, the ingrowth models on 716 and 618 observations, and the survival models on 7823 and 5244 observations, respectively, for P. sylvestris and P. nigra. The relative biases for the height models are 6.7% for P. sylvestris and 3.3% for P. nigra. The biases for the diameter growth models are zero due to the applied Snowdon correction. The biases of the ingrowth models are zero due to the applied fitting method. The relative RMSE values for the P. sylvestris and P. nigra models, respectively, are 56.4% and 48.6% for diameter growth, 24.0% and 21.7% for height, and 224.3% and 257.3% for ingrowth. 1. INTRODUCTION Pinus sylvestris L. and Pinus nigra Arn. ssp. salmannii var. pyrenaica mixtures form large forests in the Montane-Mediterranean vegetation zones of Catalonia (from 600 to 1600 m a.s.l.) [4, 32] occupying an area of 267 000 ha [12, 13]. Both species supply important products such as poles, saw logs and construction timber. The ecological (e.g. biodiversity maintenance, soil protection) and social (e.g. recreation, rural tourism, mushroom collection) functions of the pine mixtures are also significant. Most of the stands are managed using the selection system, which leads to considerable within-stand variation in tree age [11]. P. sylvestris is clearly a light demanding species, while P. nigra shows a moderate degree of shade tolerance [30], being more adaptable to irregular and multi-layered stand structures. Management planning methods currently applied in Catalonia predict the yields of stands based on yield tables and increment borings. Yield tables are static models assuming that all stands are fully stocked, pure and even-aged. They do not portray the actual or historical development of individual stands [5]. Increment borings in inventory plots are used to develop simple compartment-wise models to express diameter growth as a function of diameter. These models cannot be used in long-term simulations. Forest management planning requires growth and yield models that provide a reliable way to examine the effects of silvicultural and harvesting options, to determine the yield of each option, and to inspect the impacts of forest management on the other values of the forest [38]. Growth and yield models can be classified into two major categories: whole stand and individual tree models. Whole stand models use stand parameters such as basal area, volume, and parameters characterising the underlying diameter distribution to simulate the stand growth and yield. Individual tree models use individual trees as the basic unit for simulating tree establishment, growth and mortality; stand level values are calculated by adding the individual tree estimates together [27]. The benefit of using individual-tree models is that the stand can be illustrated much more thoroughly and several treatments simulated more easily than with stand models [29]. Individual-tree models can be distance-dependent or distanceindependent. The high cost of obtaining tree coordinates restricts the application of distance-dependent individual-tree models. The expense of such a detailed methodology is seldom warranted, making non-spatial models a more feasible alternative [38]. To date, the only empirical individual-tree growth and yield model available for the Catalan region is the non-spatial model for even-aged Scots pine stands in northeast Spain, developed by Palahí et al. [25]. Some variables such as dominant height, stand age and site index used in even-aged models are not directly applicable to uneven-aged stands [27]. The age of individual trees of an uneven-aged stand is often unknown, which means that neither stand nor tree age is a useful model predictor. An alternative to the use of these variables is to obtain site information from topographic descriptors such as elevation, slope, aspect, location descriptors (latitude), and soil type [2]. Examples of this type of models are PROGNOSIS [36, 39], designed for the Northern Rocky Mountains, PROGNAUS [22] developed for the Austrian forests, and the model developed by Schröder et al. [33] for maritime pine trees in northwestern Spain. An interesting feature of these models is that they may be applied to both uneven-aged and even-aged conditions. Another possibility to accommodate site in the model is to rely on the presence of plant species that indicate site fertility [3]. This study aims at developing a model set, which enables tree-level distance-independent simulation of the development of uneven-aged mixtures of P. sylvestris and P. nigra in Catalonia. The system consists of a diameter growth model, a static height model, an ingrowth model and a survival model for the coming 10-year period. Separate models are developed for P. sylvestris and P. nigra. The predictor variables have been restricted to site, stand and tree attributes that can be reliably obtained from stand inventories normally carried out in the region. The model set should apply to any age structure and degree of mixture (including pure stands) of the two pine species. 2. MATERIALS AND METHODS 2.1. Data The data were provided by the Spanish National Forest Inventory [6, 16–19]. This inventory consists of a systematic sample of permanent plots distributed on a square grid of 1 km, with a 10-year remeasurement interval. From the inventory plots over the whole of Catalonia, 922 plots representing all degrees of mixture (including pure stands) between P. sylvestris and P. nigra were selected (Fig. 1). The criterion for plot selection was that the occupation of one (pure stands) or two (mixed stands) of the studied species in the stands should be at least 90%. Most of the stands were naturally regenerated. The sample plots were established in 1989 and 1990. The remeasurement was carried out in 2000 and 2001. A hidden plot design was used: plot centres were marked by an iron stake buried underground; the iron stake was relocated by a metal detector. Trees were recorded by their polar coordinates and marked only temporarily during the measurements. The sampling method used circular plots in which the plot radius depended on the tree’s diameter at breast height (dbh, 1.3 m) (Tab. I). At each measurement, the following data were recorded from every sample tree: species, dbh, total height, and distance and azimuth from the plot centre. In the second measurement, a tree previously measured in the first measurement was identified as: standing, dead or thinned. Trees that entered the first dbh-class (from 7.5 to 12.4 cm) during the growth period were also recorded. The standing and dead trees resulted in 8173 diameter/ 2.2. Diameter increment modelling dbh Plot radius, m 75 ≤ dbh < 125 mm 125 ≤ dbh < 225 mm 225 ≤ dbh < 425 mm height and 8058 diameter growth observations for P. sylvestris (Tab. II), and 5721 diameter/height and 5695 diameter growth observations for P. nigra (Tab. III). There were also 721 diameter/height and 717 diameter growth observations for other species, referred to as accompanying species. Because it was not known whether a tree removed in thinning was living or dead, the thinned trees were not used as observations. At each measurement the growing stock characteristics were computed from the individual-tree measurements of the plots. A diameter growth model was prepared for both pine species. The predicted variable in the diameter growth models was the logarithmic transformation of 10-year diameter growth. This resulted in a linear relationship between the dependent and independent variables, and enabled the development of multiplicative growth models [9, 15, 21, 22, 33, 39]. Ten-year diameter growth was calculated as a difference between the two existing diameter measurements (years 1989–1990 and 2000–2001). The growth observations (10 to 12 year growth) were converted into 10-year growths by dividing the diameter increment by the time interval between the two measurements and multiplying the result by 10. The predictors were chosen from tree, stand and site characteristics as well as their transformations. All predictors had to be significant at the 0.05 level, and the residuals had to indicate a non-biased model. Due to the hierarchical structure of the data (trees are grouped into plots, and plots are grouped into provinces), the generalised least-squares (GLS) technique was applied to fit the mixed linear models. The residual variation was therefore divided a N: number of observations at tree- and stand-level; id10: 10-year diameter increment; dbh: diameter at breast height; BALsyl: competition index of P. sylvestris; BALnig+acc: competition index of P. nigra and accompanying species; BALthin: 10-year thinned competition; G: stand basal area; h: tree height; ulk: random between-plot factor; ELE: elevation; SLO: slope; LAT: latitude; CON: continentality; ING: stand ingrowth; Gsyl: stand basal area of P. sylvestris; DIN: mean dbh of ingrowth trees; P (survive): probability of a tree surviving; BALall: competition index calculated from all species. Variablea Diameter growth model (Eq. ( 1 )) id10 (cm/10 a) dbh (cm) BALsyl (m2 ha–1) BALnig+acc (m2 ha–1) BALthin (m2 ha–1) G (m2 ha–1) Diameter growth plot factor models (Eq. ( 3 )), u lk (ln (cm/10 a)) ELE (100 m) SLO (%) N into between-province, between-plot and between-tree components. The linear models were estimated using the maximum likelihood procedure of the computer software PROC MIXED in SAS/STAT [31]. The P. sylvestris (Eq. ( 1 )) and P. nigra (Eq. ( 2 )) diameter growth models were as follows: 1 ln (id10lkt)= β0 + β1 × -d---b---h---l--k--t + β2 × ln (dbhlkt) BALsyllk BALnig + acclk + β3 × ---------------------------------- +β4 × -------------------------------------ln (dbhlkt + 1) ln (dbhlkt + 1) ( 1 ) BALthinlk + β5 × ---------------------------------- + β6 × ln (Glk) + ul + ulk + elkt ln (dbhlkt + 1) 1 ln (id10lkt) = β0 + β1 × -d---b---h---l--k--t + β2 × ln (dbhlkt) BALniglk BALsyl + acclk + β3 × ---------------------------------- +β4 × ------------------------------------ ln (dbhlkt + 1) ln (dbhlkt + 1) + β5 × -----B----A----L----t--h---i--n---l--k----- + ul + ulk + elkt ( 2 ) ln (dbhlkt + 1) where id10 is future diameter growth (cm in 10 years); dbh is diameter at breast height (cm), BALsyl is the total basal area of P. sylvestris trees larger than the subject tree (m2 ha–1); BALnig + acc is the total basal area of trees that are not P. sylvestris and are larger than the subject tree (m2 ha–1); BALnig is the total basal area of P. nigra trees larger than the subject tree (m2 ha–1); BALsyl + acc is the total basal area of trees other than P. nigra and larger than the subject tree (m2 ha–1); BALthin is the total basal area of trees larger than the subject tree and thinned during the next 10-year period (m2 ha–1); and G is stand basal area (m2 ha–1). Subscripts l, k and t refer to province l, plot k, and tree t, respectively. ul, ulk and elkt are independent and identically distributed random between-province, between-plot and between-tree factors with a mean of 0 and constant variances of σp2rov, σp2l, and σt2r, respectively. These variances and the parameters βi were estimated using the GLS method. At first, all three random factors were included in the model but the between-province factor was not significant, and it was therefore excluded from the models. The random plot factors (ulk) of the models (Eqs. ( 1 ) and ( 2 )) correlated logically with the site factors. In order to include the site effects in the simulations, linear models predicting the random plot factors were developed using the ordinary least squares (OLS) technique in SPSS [35] . The models for the random plot factor of P. sylvestris (Eq. ( 3 )) and P. nigra (Eq. ( 4 )) were as follows: 2 ulk = β0 + β1 × ELElk + β2 × (ELElk) + β3 × SLOlk + elk ( 3 ) ulk = β0 + β1 × ln (ELElk) + β2 × SLOlk + β3 × CONlk the model. The non-linear height models for P. sylvestris and P. nigra were as follows: β1 hlkt = --------------------------------------------------------------------------- + ulk + elkt (1 + β2 ⁄ dbhlkt + β3 / dbhl2kt) where h is tree height (m); dbh diameter at breast height (cm); β1, β2, β3 are parameters. The random plot factors ulk were modelled as a function of site variables. The models for the random plot factor for P. sylvestris (Eq. ( 6 )) and P. nigra (Eq. ( 7 )) were developed using the ordinary least squares (OLS) technique in SPSS [35]: ulk = β0 + β1 × ELElk + β2 × LATlk + β3 × CONlk 2 + β4 × ln (CONlk) + β5 × (CONlk) + elk ulk = β0 + β1 × ELElk + β2 × LATlk + β3 × (CONlk) 2 1 + β4 × ---------------- + elk CONlk + β4 × LATlk + elk where ulk is plot factor predicted by equations ( 1 ) or ( 2 ); ELE is elevation (100 m); SLO is slope (%); CON is continentality (linear distance to the Mediterranean Sea, km); LAT is latitude (y UTM coordinate, 100 km). In simulations, the random plot factor (ulk in Eqs. ( 1 ) or ( 2 )) may be replaced by its prediction (Eqs. ( 3 ) or ( 4 )). Other site characteristics and their transformations adopted logical signs, namely aspect, soil texture, and humus, but were not significant. Another version of the plot factor models was prepared using the presence of certain plant species in the stand as dummy variables (referred to as species dummies), in addition to variables listed in equations ( 3 ) and ( 4 ). To convert the logarithmic predictions of equations ( 1 ) and ( 2 ) to the arithmetic scale, a multiplicative correction factor suggested by Baskerville [1] was tested (exp(s2/2)), where s2 is the total residual variance of the logarithmic regression). However, it resulted in biased back-transformed predictions. Therefore, an empirical ratio estimator for bias correction in logarithmic regression was applied to equations ( 1 ) and ( 2 ). As suggested by Snowdon [34], the proportional bias in logarithmic regression was estimated from the ratio of the mean diameter growth id10 and the mean of the back-transformed predicted values from the regression exp [ ln iˆd 10] . The ratio estimator was therefore id10 ⁄ exp [ ln iˆd 10] . 2.3. Height modelling Analysis of the height data revealed that there were obvious and large errors in the height measurements of the first measurement occasion. Therefore, height growth models could not be estimated. Consequently, static height models using the second measurement were developed. Models that enable the estimation of total tree heights when only tree diameters and site characteristics are measured (as is the case in forest inventory) were estimated. Elfving and Kiviste [8] proposed 13 functions having a zero point, being monotonously increasing and having one inflexion point, for approximation of the relationship between stand age and height. These functions were tested as the height model, but dbh was used instead of age as the predictor. A total of 10 two- and three-parameter functions were tested. The models developed by Hossfeld [28] and Verhulst [14] gave the best fit. Out of this these, Hossfeld model (Eq. ( 5 )) was selected because it has been used earlier in Spain [24, 26]. The non-linear height models were estimated using the non-linear mixed procedure (NLMIXED) in SAS/STAT [31]. In the procedure, it is possible to include only two random factors in the model. Because the random between-plot factor was more significant than the random between-province factor, the plot factor was included in ( 4 ) where ulk is random plot factor of the related height model. Other site characteristics and their transformations such as aspect, slope and soil texture were not significant in the final version of the models. Another version of the models was prepared using species dummies as additional predictors. 2.4. Ingrowth modelling A linear model predicting the number of trees per hectare entering the first dbh-class (from 7.5 to 12.4 cm) during a 10-year growth period was prepared for each species. The predictors were chosen from stand and site characteristics and their transformations. Mixed linear models were estimated first, but the random between-province factor was not significant. Thus, ingrowth models for P. sylvestris (Eq. ( 8 )) and P. nigra (Eq. ( 9 )) were estimated using the ordinary least squares (OLS) method in SPSS [35]: 1 Gsyl INGlk = β0 + β1 × Glk + β2 × ------- + β3 × -------------l-k- + elk Glk Glk Gnig CON INGlk = β0 + β1 × Glk + β2 × -----G----l--k---l--k + β3 × -E----L----E----l-l-kk- + elk where ING is ingrowth (number of trees ha–1) at the end of a 10-year growth period; Gsyl and Gnig are stand basal area of P. sylvestris and P. nigra, respectively (m2 ha–1). The mean dbh of the ingrowth trees of P. sylvestris (Eq. ( 10 )) and P. nigra (Eq. ( 11 )) was modelled as well. The models were estimated using the ordinary least squares (OLS) method: DINlk = β0 + β1 × Glk + β2 × ELElk + elk DINlk = β0 + β1 × Glk + elk where DIN is the mean dhh of ingrowth trees (cm) at the end of a 10year growth period. Another version of the models using species dummies as predictors for the number and mean dbh of ingrowth was also evaluated. 2.5. Survival modelling When analysing the data, two types of mortality were identified: density-independent mortality and density-dependent. The densityindependent tree-level survival rate for a 10-year period was estimated at 0.962 overall. All mortality of plots having basal area values at the second measurement lower than 1 m2 ha–1 or lower than 90% of the stand basal area at the first measurement were considered as density-independent (usually caused by fire). A model for the density-dependent probability of a tree to survive for the next 10-year growth period was estimated from the remaining sample plots. The following survival models for P. sylvestris ( 5 ) ( 6 ) ( 7 ) ( 8 ) ( 9 ) ( 10 ) ( 11 ) 1 P(survive)lkt = ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- + elkt 1 + exp –β0 + β1 × ----------B----A----L---l--k--t--------- + β2 × hlkt + β3 × ELElk + β4 × CONlk ln (dbhlkt + 1) 1 P(survive)lkt = ------------------------------------------------------------------------------------------------------------------------------------------------------------- + elkt 1 + exp –β0 + β1 × ----------B----A----L---l--k--t--------- + β2 × Glk + β3 × ELElk ln (dbhlkt + 1) ( 12 ) ( 13 ) (Eq. ( 12 )) and P. nigra (Eq. ( 13 )) were estimated using the Binary Logistic procedure in SPSS [35]. See equations (12) and (13) above where P(survive) is the probability of a tree surviving for the next 10year growth period. Another version of the models was developed using the presence of particular plant species as a site fertility indicator. 2.6. Model evaluation 2.6.1. Fitting statistics The models were evaluated quantitatively by examining the magnitude and distribution of residuals for all possible combinations of variables included in the model. The aim was to detect any obvious dependencies or patterns that indicate systematic discrepancies. To determine the accuracy of model predictions, the bias and precision of the models were calculated [10, 21, 25, 38]. The absolute and relative biases and the root mean square error (RMSE) were calculated as follows: ∑ (yi – yˆi) bias = --------------------------n ∑ (yi – yˆ i) / n bias% = 100 × ---------------------------------- ∑ yˆi / n RMSE = ∑ (yi – yˆi)2 ---------------------------- n – 1 RMSE% = 100 × R2 = 1 – -∑--------(--y---i---–-----yˆ---i-)--2∑ (yi – y)2 ∑ (yi – yˆi)2 ⁄ (n – 1) ----------------------------------------------------∑ yˆi / n ( 14 ) ( 15 ) ( 16 ) ( 17 ) ( 18 ) where n is the number of observations; and yi and yˆ i are observed and predicted values, respectively. In the models that included a random plot factor, the predicted value (yˆ i) was calculated using a model prediction of the plot factor. 2.6.2. Simulations In addition, the models were further evaluated by graphical comparisons between measured and simulated stand development. The simulated 10-year change in stand basal area of the inventory plots was compared to the measured change. The dynamics of accompanying species, present in several plots, was simulated using equations shown in the Appendix. The simulation of a 10-year time step consisted of the following steps: 1. For each tree, add the 10-year diameter increment (Eqs. ( 1 ) and ( 2 )) using the predicted plot factor (Eqs. ( 3 ) and ( 4 )) to the diameter. 2. Multiply the frequency of each tree (number of trees per hectare that a tree represents) by the density-dependent 10-year survival probability. The density-dependent probability is provided by equations ( 12 ) and ( 13 ). 3. Calculate the number of trees per hectare (Eqs. ( 8 ) or ( 9 )) that enter the first dbh-class and the mean dbh of ingrowth (Eqs. ( 10 ) or ( 11 )) at the end of a 10-year growth period. 4. Calculate tree heights using equation ( 5 ), and the predicted plot factor provided by equations ( 6 ) or ( 7 ). In addition, the development of two plots – one representing a mixed P. sylvestris and P. nigra stand and another representing a pure stand of P. sylvestris – was simulated at different elevations to evaluate the model set in long-term simulation. 3. RESULTS 3.1. Diameter growth models Parameter estimates of the diameter growth models (Eqs. ( 1 ) and ( 2 )) were logical and significant at the 0.001 level (Tab. IV). Parameter estimates of the plot factor models were significant at the 0.05 level. The R2 values were 0.13 and 0.14 for the P. sylvestris and P. nigra diameter growth models, respectively. The R2 value of the random plot factor model was 0.06 for P. sylvestris and 0.10 for P. nigra, showing that only a small part of the variation in plot factor was explained by site characteristics. The explained variation was higher when species dummies were used, resulting in R2 values of 0.11 for P. sylvestris and 0.18 for P. nigra. The R2 values of predictions using both the diameter growth and plot factor models (Eq. ( 18 )) were 0.16 for P. sylvestris and 0.18 for P. nigra. When using species dummies in the plot factor models, these values were 0.18 for P. sylvestris and 0.21 for P. nigra. The shape of the relationship between dbh and diameter growth is typical of tree growth processes ([39], Fig. 2). Diameter increment of dominant trees (BALx = 0) increases to a maximum at dbh of 17 cm and then slowly decreases, approaching zero asymptotically as the tree matures (Eqs. ( 1 ) and ( 2 )). Increasing competition (G, BALsyl and BALnig + acc in Eq. ( 1 ); BALnig and BALsyl + acc in Eq. ( 2 )) decreases the diameter growth. The models indicate that P. nigra causes more competition because the coefficients of competition calculated from P. nigra trees (β4 in Eq. ( 1 ) and β3 in Eq. ( 2 )) always had higher absolute values than BAL computed from P. sylvestris (β3 in Eq. ( 1 ) and β4 in Eq. ( 2 )). The thinned competition (BALthin) had a positive effect on diameter growth (Eqs. ( 1 ) and ( 2 )) (Fig. 3). This variable improved the fit and logical behavior of the other predictors in the models, although the variable is seldom used when the models are applied in simulation (i.e. this variable is given a zero value). Increasing slope decreased the plot factor and consequently the diameter growth of all trees on a plot (Eqs. ( 3 ) and ( 4 )). According to the models, elevation affects differently the two studied species: higher growth rates of P. sylvestris are observed at extreme elevations (Fig. 4), while a S.E. of estimates are given in parenthesis. b ROS is Rosa spp., JPH is Juniperus phoenicea, ROM is Rosmarinus officinalis, CRA is Crataegus sp., JUN is Juniperus communis. increasing elevation increases the growth of a P. nigra tree. The signs of coefficients of the plot factor model of P. nigra were logical, bearing in mind the climatic models (predicting mean extreme temperatures and precipitation) developed by Ninyerola et al. [23] for the Catalan region: increasing continentality decreases the growth of a tree, and the more northern the latitude the higher is the stand growth (Fig. 5). The ratio estimators for bias correction in the fixed part of the P. sylvestris and P. nigra diameter growth models (Eqs. ( 1 ) and ( 2 )) were 2.6324/2.1288 = 1.2365 and 2.7981/2.5352 = 1.1037, respectively. The ratio estimators for bias correction using both the fixed part and the predicted plot factors (Eqs. ( 1 ), ( 2 ), ( 3 ) and ( 4 )) were 2.6324/2.1389 = 1.2307 for P. sylvestris and 2.7981/2.5639 = 1.0914 for P. nigra. When using species dummies in the plot factor models, the ratio estimators were 2.6324/2.1453 = 1.2270 for P. sylvestris and 2.7981/2.4847 = 1.1261 for P. nigra. The bias of the growth models, when the fixed model part and the plot factor models without species dummies were used, showed no trends when displayed as a function of predictors or predicted growth in Figures 6 and 7. The residuals of the diameter growth and height models are correlated within each plot – part of the residual variation is explained by random between-plot factor, but only a small part of the between-plot variation is explained by plot factor model. This should be taken into account when analyzing Figures 6 and 7. The absolute and relative biases for P. sylvestris and P. nigra diameter growth models were zero due to the ratio estimator used for bias correction. The relative RMSE values were 56.4% and 48.6% for the P. sylvestris and P. nigra models, respectively (Tab. V and VI). 3.2. Height models The estimated height models describe tree height as a function of diameter at breast height (Eq. ( 5 )). According to the models for the random plot factor (Eqs. ( 6 ) and ( 7 )), the effect of site characteristics on tree height of both pine species is very similar: increasing elevation decreases the height of a tree; continentality first increases the height (up to around 80 km) and then decreases the height; and latitude increases the height. The use of species dummies resulted in a clear improvement of the plot factor models. Parameter estimates of the height models and plot factor models were logical and significant at the 0.05 level (Tab. VII). The R2 values were 0.30 for the P. sylvestris height model, 0.41 for the P. nigra height model, 0.15 for the P. sylvestris plot factor model without species dummies, 0.18 for the P. nigra plot factor model without species dummies, 0.32 for the P. sylvestris plot factor model with species dummies, and 0.35 for the P. nigra plot factor model with species dummies. The R2 values, when adding the predicted plot factor to the fixed part of the height model, were 0.33 (0.43 using species dummies) for P. sylvestris and 0.47 (0.55 using species dummies) for P. nigra. The relative biases were 6.7% and 3.3% and the relative RMSE were 24.0% and 21.7% for the P. sylvestris and P. nigra height models, respectively (Tabs. V and VI). There were no obvious trends in bias for the height models, but the residuals had a slightly heterogeneous variance as a function of predicted height. CRA ACR FAG THI JUN σ2pl σ2tr RMSE R2 Parameter β0 β1 β2 β3 σp2l R2 (Tab. VIII). The R2 values were 0.11 for the P. sylvestris ingrowth model, 0.11 for the P. nigra ingrowth model, 0.12 for the P. sylvestris mean dbh of ingrowth model, and 0.05 for the P. nigra mean dbh of ingrowth model. The developed models a S.E. of estimates are given in parenthesis. b CRA is Crataegus sp., ACR is Acer sp., FAG is Fagus sylvatica, THI is Thimus ssp., JUN is Juniperus communis. a S.E. of estimates are given in parenthesis. 3.3. Ingrowth models Parameter estimates of the models for the number and mean dbh of ingrowth were logical and significant at the 0.05 level for the number and mean dbh of ingrowth use stand basal area and site characteristics as independent variables: increasing stand basal area increases the amount of ingrowth for P. sylvestris (Eq. ( 8 )) up to 8 m2 ha–1, and then decreases the amount of ingrowth; increasing stand basal area decreases the amount of ingrowth for P. nigra (Eq. ( 9 )) and the mean dbh of ingrowth for both pine species (Eqs. ( 10 ) and ( 11 )); increasing values of the ratio of the subject species’ basal area to the total stand basal area increases ingrowth (Eqs. ( 8 ) and ( 9 )); increasing ratio between continentality and elevation increases P. nigra ingrowth; and increasing stand elevation increases the mean dbh of P. sylvestris ingrowth (Eq. ( 10 )). The main difference between both species relies on higher ingrowth for P. nigra, at low and continental sites (β3 in Eq. ( 9 )). The use of species dummies did not bring about a significant improvement in the models. The absolute and relative bias for the P. sylvestris and P. nigra ingrowth and mean dbh of ingrowth models were zero. The relative RMSE value was 224.3% for P. sylvestris ingrowth, 257.3% for P. nigra ingrowth, 10.1% for P. sylvestris mean dbh of ingrowth, and 10.2% for P. nigra mean dbh of ingrowth. There were no obvious trends in the bias for the ingrowth models, but the residuals had a slightly heterogeneous variance as a function of predicted ingrowth, P. sylvestris basal area, P. nigra basal area and elevation. The graphs of bias as a function of predicted ingrowth showed that the models may predict negative ingrowth (e.g. with very high stand basal area). In simulations, negative ingrowth predictions should be replaced by zero. β0 β1 β2 β3 β4 CIT CAV JUN COA ACR BUX χ2-value β0 β1 β2 β3 χ2-value β0 β1 β2 β3 PRU AUU χ2-value P. sylvestris survival model with sp. dummies (Eq. ( 12 )) P. nigra survival model without sp. dummies (Eq. ( 13 )) P. nigra survival model with sp. dummies (Eq. ( 13 )) a CIT is Cistus spp., CAV is Calluna vulgaris, JUN is Juniperus communis, COA is Corylus avellana, ACR is Acer sp., BUX is Buxus sempervirens, PRU is Prunus sp., AUU is Arctostaphylos uva-ursi. 3.4. Survival models The probability of a tree surviving for the next 10 years was best explained by the ratio of the basal area of trees larger than the subject tree to the subject tree’s dbh (BAL/ln(dbh+1)) (Eqs. ( 12 ) and ( 13 )), tree height (Eq. ( 12 )), continentality (Eq. ( 12 )), stand basal area (Eq. ( 13 )), and elevation (Eqs. ( 12 ) and ( 13 )). The Wald tests showed that the parameter estimates of equations ( 12 ) and ( 13 ) are significant (P < 0.05) (Tab. IX). By analysing equations ( 12 ) and ( 13 ) it can be deduced that: ( 1 ) the greater the ratio of basal area of trees larger than the subject tree (competition index) to tree dbh, the smaller the survival probability; ( 2 ) the taller the tree and the higher the elevation, the greater the probability of a P. sylvestris tree to survive; ( 3 ) the higher the elevation, the smaller the survival probability of a P. nigra tree; ( 4 ) the denser the stand, the greater the survival probability of a P. nigra tree; and ( 5 ) the greater the continentality, the smaller the probability of a P. sylvestris tree surviving. The main differences between the two species rely on survival rate for short trees (β2 in Eq. ( 12 )) and dense stands (β2 in Eq. ( 13 )). The odds ratios of the covariates showed that BAL/ ln(dbh+1) (Eqs. ( 12 ) and ( 13 )) has the strongest relative effect on the probability of a tree surviving. The use of species dummies gave a significant improvement of the survival models [7]. 3.5. Simulation results Figure 8 shows long-term simulations without any treatments of two plots, one representing a mixed P. sylvestris and P. nigra stand, and the other representing a pure stand of P. sylvestris. The simulations were carried out for four different elevations (200, 500, 1500, 1900 m a.s.l.). Site characteristics were the same for both plots. The long-term simulations include the density-independent survival of 0.962 per 10 years. The predicted plot factors (without using species dummies) have been added to the predictions and the diameter increment models have been corrected for logarithmic transformation. In simulations, P. nigra occupies the site at lower elevations in both plots, and even tends to dominate at 1500 a.s.l. in the mixed plot. At higher elevations (from 1500 m a.s.l. upwards in P. sylvestris, and at 1900 m a.s.l. in the mixed plot) P. sylvestris occupies the stands. The simulation results correspond to the observed species distribution. Maximum volume is higher at low elevations. However, the volume stops decreasing at higher elevations (1500 m a.s.l.) and even increases to some extent at the highest elevation (1900 m a.s.l.). Figure 9 shows the measured and predicted 10-year changes of stand basal area for those non-thinned plots that were used to develop the density-dependent survival models. In the simulations the density-independent 0.962 survival per 10 years was not applied, but the diameter increment models were corrected for logarithmic transformation. In Figure 9A just the fixed part of the diameter increment model was used; in Figure 9B the random plot factors predicted without species dummies were added to the fixed part of the models; in Figure 9C the random plot factors predicted using species dummies were added to the fixed part of the models; and in Figure 9D the fixed part of the models and the true random plot factors were used. There was a small improvement in the predictions when the predicted plot factors were used. The predictions were slightly better when the predictions were based on models that use species dummies. The improvement was more significant when the true plot factors were added to the fixed part of the model. The model set underpredicted the 10-year change in stand basal area of plots having exceptionally high growth. A careful inspection of the plot data showed that plots with the highest underprediction were young fastgrowing, rather even-aged stands. The residuals of the simulated 10-year change in stand basal area were also plotted against site characteristics (results not shown). No systematic trends were found, showing that there were no model failures related to these variables in the models. Nonetheless, the residuals were positively biased due to the inability of the model set to predict high enough growth in young fast-growing stands. 4. DISCUSSION This study presents individual-tree models for uneven-aged mixtures of P. sylvestris and P. nigra in Catalonia, based on permanent sample plots measured two times in all sites represented by the Spanish National Forest Inventory. This sample provides an outstanding database in terms of size (14 470 diameter growth observations) and forest conditions. However, it should be taken into account that the sampling methodology was not specifically designed to develop growth and yield models. A disadvantage of this data is that diameter growth is determined as a difference of two diameters. The breast height diameter may not have been measured at exactly the same height, and the direction of the diameter measurement may have been different on different measurement occasions. This results in greater errors than measuring radial increment directly from increment borings. This is reflected in the value of the coefficients of determination, 0.18 for P. sylvestris and 0.21 for P. nigra, when the species dummies are used in the plot factor models. The low R2 of this study agrees with the results obtained by Monserud and Sterba [22], using similar data from the Austrian National Forest Inventory. Nevertheless, assuming that the measurement errors were random, the large sample should compensate for this disadvantage. The variable-radius circular plot sampling method, used to collect the modelling data, selected trees with unequal probability (Tab. I). This method results in ingrowth to dbh classes larger than the smallest class, 7.5–12.4 cm. The ingrowth models prepared in this study predict only the number of trees that enter the smallest dbh class during a 10-year period. However, despite its specific features, the variable radius circular plot is often the only feasible method to sample irregular stand structures efficiently. This sampling method gives a good representation of large trees, which is usually a benefit from both inventory and modelling stand points. Sampling methods often have limitations to represent spatial variability in stands [37]. Hence, competition predictors used in the models might have sampling error associated with them, which will create bias when using the models in simulations. The site effects were modelled from site characteristics (e.g. elevation), avoiding the use of site index and stand age. More accurate description of site fertility was obtained when species dummies were used, although it should be remembered that the presence of species in stands can be affected by management, forest fires or grazing. The between-plot factors in the diameter growth and height models correlated logically with site variables. Nevertheless, site characteristics adopted illogical signs when included directly in the fixed part of the models. For that reason, site characteristics were used to develop models that predict the plot factors of the diameter increment and height models. Because it was not known whether a tree removed in thinning was living or dead, the thinned trees were not used as observations. If dead thinned trees had been recorded differently from live trees, the survival and growth models might have been better for trees facing much competition. Height growth models were not developed because there were obvious and large errors in the height measurement of the first measurement occasion. The developed height prediction models provide typical values for a given dbh-class and can be used for growth simulation. The number and mean dbh of ingrowth can be predicted using stand level models. A density-independent 10-year tree survival rate of 0.962 was estimated from the study material. However, this rate should not be automatically used in simulation because the rate may depend on location, fire management, etc. In most cases it might be better to use only the density-dependent survival models. Models for the occurrence of fire may improve the prediction of density-independent mortality. Long-term simulations with the model set suggest that P. nigra occupies the site at lower elevations if no treatments are applied, which is reasonable if we take into account that this species shows higher shade tolerance than P. sylvestris. Furthermore, succession in temperate forests appears to be driven by differences in light availability and shade tolerance; however, water limitation is also important for the distribution of forest species in Mediterranean plant communities [40]. The ingrowth and survival models do not have a strong impact on short-term simulations, but in long-term simulations they are very important. The system of models shows some weaknesses with very large trees (over 50 cm of dbh), because the data included very few large trees. The growth of very large trees may be overpredicted and the tree survival models may underpredict the mortality of very large trees. Other authors [15, 21, 33, 39] used squared dbh when modelling the effect of tree size, in order to get the growth predictions to approach zero for large dbh values. However, dbh2 was not a significant predictor in our data. Hence, long-term simulations may produce trees older than 200–300 years, which is the observed range of maximum tree age (244 years for P. sylvestris and 215 years for P. nigra), provided by the Ecological and Forest Inventory of Catalonia [12, 13]. An age-dependent survival model could be considered in future studies. However, the weaknesses of the models with very large trees do not have any practical importance in short- and medium-term simulations because very large trees are rarely found in these types of stands. The models provide correct average predictions (Fig. 9), but they account for only a small part of the site specific growth variation among stands. The highest underpredictions are related to young fast-growing, rather even-aged stands. The lack of age and site index reduces the ability of the models to accurately predict growth in even-aged stands. Nevertheless, the unexplained site-specific variation could be accounted for by using model calibration [15, 20, 22, 39]. Figure 9D, where the true plot factors have been used, gives an idea about the performance of a calibrated model. The potential for application of the models is wide. They can be applied using the data normally available from stand inventories in the region, and they can be adapted to all age structures and degrees of mixture, including pure stands. However, as mentioned above, for even-aged stands, models that use site index as a predictor should be considered as the first option. This study being the first, known by the authors, on individual-tree growth models for uneven-aged mixed stands in Spain is a starting point for further research on the topic. Future studies could evaluate the effect of using other potential predictors in the models, such as tree age, crown width or length (characterising the variation in the vigour of trees of similar size) [22, 33, 39], or additional locally measured site variables, such as soil depth [22, 33]. Future studies should also focus on evaluating the use of past increment as a model predictor or for calibrating the models for a specific stand. The models presented in this study can be used to optimise the stand management and to evaluate alternative management regimes for unevenaged stands of P. sylvestris and P. nigra. Acknowledgements: Financial support for this project was provided by the Forest Technology Centre of Catalonia (Solsona, Spain). We are grateful to Jose Antonio Villanueva, head of Spanish National Forest Inventory, for making the Forest Inventory data available. We thank him and the staff of the “Inventario Forestal Nacional” for their cooperation and assistance. We also thank Carlos Gracia and Jordi Vayreda from the Ecological and Forest Inventory of Catalonia, for allowing us to use their data and for helping us to assess different aspects of it. We thank Associate Editor Prof. David E. Hibbs and the anonymous reviewers for their valuable comments. We thank Mr Tim Green for the linguistic revision of the manuscript. APPENDIX In addition to P. sylvestris and P. nigra, other tree species where also present in several stands. The accompanying species were divided into two major groups: conifers and hardwoods. A set of models was developed for each group of accompanying species using the ordinary least squares (OLS) method in SPSS (SPSS Inc., 1999), in order to include them in simulations. Conifers: ln (id10lkt) = 7.680 – 32.177 ⁄ dbhlkt – 1.495 × ln (dbhlkt) – 0.020 × BALacclk – 0.052 × ELElk + elkt 36.755 – 2.407 × ( cos (ASPlk) × ln (ELElk)) hlkt = --------------------------------------------------------------------------------------------------------------1 + (43.405 / dbhlkt) – 0.046 × CONlk + 2.852 × sin (ASPlk) + ------------------------1-----+----(---4---3---.-4---0---5-----/----d---b---h----l--k---t--)------------------------ + elkt SVlk = 0.9939 INGlk = 5.28 Hardwoods: SV = 0.9739 lk ING = 32.77 lk DIN = 6.27 lk – 0.004 × CON – 0.177 × sin (ASP ) + e lk lk lkt 13.644 – 0.231 × ELE + 0.023 × CON lk lk h = -------------------------------------------------------------------------------------------------------- + e lkt lkt 2  1 + 130.795 / dbh  lkt where id10 is future diameter growth (cm per 10 years); dbh is diameter at breast height (cm); G is stand basal area (m2 ha–1); BALacc is the total basal area of accompanying tree species larger than the subject tree (m2 ha–1); ELE is elevation (100 m); CON is continentality (linear distance to the Mediterranean sea, km); ASP is aspect (rad); SV is survival rate for a 10-year time step; ING is ingrowth (number of trees ha–1) at the end of a 10-year time step; DIN is the mean dhh of ingrowth trees (cm) at the end of a 10-year step. The ratio estimators for bias correction in logarithmic diameter growth models were 1.3220 and 1.2386 for conifers and hardwoods, respectively. [31] SAS Institute Inc., SAS/STAT User’s guide, version 8, Cary, NC: SAS Institute Inc., 1999, 3884 p. [35] SPSS Inc., SPSS Base system syntax reference Guide. Release 9.0, 1999. [36] Stage A.R., Prognosis model for stand development, USDA For. Ser. Res. Pap. INT-137, 1973, p. 32. [37] Stage A.R., Wykoff W.R., Adapting distance-independent forest growth models to represent spatial variability: effects of sampling design on model coefficients, For. Sci. 44 (1998) 224–238. [39] Wykoff R.W., A basal area increment model for individual conifers in the northern Rocky Mountains, For. Sci. 36 (1991) 1077–1104. [1] Baskerville G.L. , Use of logarithmic regression in the estimation of plant biomass , Can. J. For. Res . 2 ( 1972 ) 49 - 53 . [2] Bravo F. , Montero G. , Site index estimation in Scots pine (Pinus sylvestris L.) stands in the High Ebro Basin (northern Spain) using soil attributes , For . 74 ( 2001 ) 395 - 406 . [3] Cajander A.K. , The theory of forest types , Acta For. Fenn . 29 ( 1926 ) pp. 108 . [4] Conesa Mor J.A. , Tipologia de la vegetació: anàlisi i caracterització , Universitat de Lleida, 1997 . [5] Davis L.S. , Johnson K.N. , Forest management , McGraw-Hill , New York, 1987 . [6] DGCN , Tercer Inventario Forestal Nacional (1997-2006) Galicia : A Coruña, Ministerio de Medio Ambiente, Madrid, 2001 . [7] Eid T. , Tuhus E. , Models for individual tree mortality in Norway, For . Ecol. Manage. 154 ( 2001 ) 69 - 84 . [8] Elfving B. , Kiviste A. , Construction of site index equations for Pinus sylvestris L. using permanent plot data in Sweden, For . Ecol. Manage. 98 ( 1997 ) 125 - 134 . [9] Flewelling J.W. , Pienaar L.V. , Multiplicative Regression with Lognormal Errors, For . Sci. 27 ( 1981 ) 281 - 289 . [10] Gadow K. , von Hui G. , Modeling forest development , Faculty of Forest and Wodland Ecology , University of Göttingen, 1998 . [11] Gonzalez J.M. , Arrufat D. , Meya D. , Modelos de gestión selvícola para masas irregulares de pino laricio en el Prepirineo Catalan , Revista Forestal Española 16 ( 1997 ) 14 - 20 . [12] Gracia C. , Abril M. , Barrantes O. , Burriel J.A. , Ibàñez J.J. , Serrano M.M. , Vayreda J. , Inventari Ecològic i Forestal de Catalunya: Métodes, Departament d'Agricultura, Ramaderia i Pesca , Generalitat de Catalunya, Barcelona, 1992 . [13] Gracia C. , Burriel J.A. , Ibàñez J.J. , Mata T. , Vayreda J. , Inventari Ecològic i Forestal de Catalunya: Regió Forestal IV , CREAF, Bellaterra, 2000 . [14] Grimm H. , On growth models and analysis of growth curves in microbiology, Biom . J. 19 ( 1977 ) 529 - 534 . [15] Hökkä H. , Groot A. , An individual-tree basal area growth model for black spruce in second-growth peatland stands , Can. J. For. Res . 29 ( 1999 ) 621 - 629 . [16] ICONA , Segundo Inventario Forestal Nacional (1986-1995) Cataluña : Barcelona, MAPA , Madrid, 1993 . [17] ICONA , Segundo Inventario Forestal Nacional (1986-1995) Cataluña : Girona, MAPA , Madrid, 1993 . [18] ICONA , Segundo Inventario Forestal Nacional (1986-1995) Cataluña : Lleida, MAPA , Madrid, 1993 . [19] ICONA , Segundo Inventario Forestal Nacional (1986-1995) Cataluña : Tarragona, MAPA , Madrid, 1993 . [20] Lappi J. , Bailey R.L. , A height prediction model with random stand and tree parameters: an alternative to traditional site index, For . Sci. 34 ( 1988 ) 907 - 927 . [21] Mabvurira D. , Miina J. , Individual-tree growth and mortality models for Eucalyptus grandis (Hill) Maiden plantations in Zimbabwe, For . Ecol. Manage. 161 ( 2002 ) 231 - 245 . [22] Monserud R.A. , Sterba H. , A basal area increment model for individual trees growing in even- and uneven-aged forest stands in Austria, For . Ecol. Manage. 80 ( 1996 ) 57 - 80 . [23] Ninyerola M. , Pons X. , Roure J.M., A methodological approach of climatological modelling of air temperature and precipitation through GIS techniques , Int. J. Climatol . 20 ( 2000 ) 1823 - 1841 . [24] Palahí M. , Tomé M. , Pukkala T. , Trasobares A. , Montero G. , Site index model for Pinus sylvestris in north-east Spain, For . Ecol. Manage. 187 ( 2004 ) 35 - 47 . [25] Palahí M. , Pukkala T. , Miina J. , Montero G. , Individual-tree growth and mortality models for Scots pine (Pinus sylvestris L.) in northeast Spain , Ann. For. Sci. 60 ( 2003 ) 1 - 10 . [26] Pita P.A. , La calidad de la estación en las masas de Pinus sylvestris de la Península Ibérica, Anales del Instituto Forestal de Investigaciones y experiencias 9 ( 1964 ) 5 - 28 . [27] Peng C. , Growth and yield models for uneven-aged stands: past, present and future , For. Ecol. Manage . 132 ( 2000 ) 259 - 279 . [40] Zavala M.A. , Constraints and trade-offs in Mediterranean plant communities: The case of Holm Oak-Aleppo Pine Forests, Bot . Rev. 66 ( 2000 ) 119 - 149 .


This is a preview of a remote PDF: http://www.afs-journal.org/articles/forest/pdf/2004/01/F4102.pdf

Antoni Trasobares, Timo Pukkala, Jari Miina. Growth and yield model for uneven-aged mixtures of Pinus sylvestris L. and Pinus nigra Arn. in Catalonia, north-east Spain, Annals of Forest Science, 9-24, DOI: doi:10.1051/forest:2003080