Development of a system to predict the evolution of individual tree mature cork caliber over time

Forest Ecology and Management 260 (2010) 1303–1314

Contents lists available at ScienceDirect

Forest Ecology and Management journal homepage: www.elsevier.com/locate/foreco

Development of a system to predict the evolution of individual tree mature cork caliber over time Alice M. Almeida ∗ , José Tomé, Margarida Tomé Centro de Estudos Florestais, Instituto Superior de Agronomia, Universidade Técnica de Lisboa, Tapada da Ajuda, 1349 – 017 Lisboa, Portugal

a r t i c l e

i n f o

Article history: Received 4 April 2010 Received in revised form 6 July 2010 Accepted 10 July 2010 Keywords: Cork caliber Cork thickness of complete rings Cork growth models Cork growth index

a b s t r a c t The development of a model for the prediction of the evolution of individual tree cork caliber over time, from a measurement taken at a certain point in time, was the main objective of this work. The model includes three sub-models: a model to predict the thickness of complete rings from cork caliber at tc years; a cork growth model (for complete rings) and a model to predict cork caliber at age tc from the corresponding thickness of (tc − 1) complete rings. The algebraic difference approach (ADA) as well as the generalized algebraic difference approach (GADA) were used in modeling cork growth. Several base models with one or two site-tree-specific parameters were fitted to the data using the dummy variable approach. The selection of the cork growth model was based on several criteria: fitting ability, prediction performance evaluated through the PRESS residuals and behaviour screened with available knowledge on the cork growth process. The ADA model derived from the log-logistic function with the asymptote as free parameter was selected. The models developed to predict cork caliber and the corresponding thickness of complete rings were based in the linear relationship between the two variables. The two models were simultaneously fitted using two stage least squares approach. The predicted thickness of cork complete rings in a 9-year old cork is proposed as a cork growth index. The distribution of this index can be used to characterize the potential of a site for cork production. © 2010 Elsevier B.V. All rights reserved.

1. Introduction Cork oak (Quercus suber L.) is an evergreen oak known for its ability to produce a thick and continuous layer of suberised cells, the cork. The cork industry transforms this raw material in different products, but the economic feasibility of the whole sector is dependent on the production of stoppers of natural cork for use in the bottling of wines. The suitability of the raw material for this product establishes its commercial value and therefore the objectives of forest management. The thickness of the cork plank, or cork caliber according to the industry terminology, is the most important variable when analyzing the raw material suitability for the production of stoppers. After extraction from the tree, the cork planks undergo a postharvest preparation for further industrial processing consisting of an immersion in water at approximately boiling temperature during 1 h. The objective of this operation is to flatten the raw planks, curved according to the stem shape, and to soften the cork tissue for an easier subsequent cutting. With water boiling the cork expands and the most important practical consequence is that the raw cork

∗ Corresponding author. Tel.: +351 21 363 33 56; fax: +351 21 364 50 00. E-mail addresses: [email protected] (A.M. Almeida), [email protected] (J. Tomé), [email protected] (M. Tomé). 0378-1127/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.foreco.2010.07.017

planks increase in thickness, on average 12% (Pereira and Tomé, 2004). Cork thickness can be measured either before or after boiling but the last measurement is the reference for industry. Additionally, cork thickness after boiling best represents the real cork thickness, since it is measured after the internal tensions, caused by the cellular corrugation during cork growth, occurring between the wood and the external cork layers, have been relieved during the water boiling of cork (Pereira, 2007). This effect is particularly important in the radial direction along which cork thickness is measured. In this paper, unless otherwise indicated, the term cork thickness refers to the value measured after boiling. Cork thickness depends on environmental conditions and tree genetics that determine cork growth rate, but also on the duration of the cork production cycle, defined as the number of years between two cork extractions. The longer the period between cork extractions, the thicker the cork layer that is produced. The duration of the cycle is dictated by the need to obtain a certain thickness in the cork plank for its subsequent industrial utilization. The total cork thickness that allows the production of stoppers corresponds to a minimum of 27 mm while corks with a total cork thickness greater than 40 mm, having a greater percentage of waste (Adrados et al., 2000), must be avoided. According to Portuguese legislation, the interval between two consecutive cork extractions on the same tree must be equal or greater than 9 years. Although the majority of the cork oak stands

1304

A.M. Almeida et al. / Forest Ecology and Management 260 (2010) 1303–1314

Fig. 1. Schematic representation of a cork sample with 9 years, showing the 8 complete cork rings and the 2 half rings (adapted from Natividade, 1940).

are debarked with this periodicity, the question whether it is better to reduce or extend the interval between cork extractions according to the site and/or the market requirements arises very often. It is therefore important to have models to estimate the cork growth to help to assess the most favourable cork age for extraction in each stand in order to get corks that maximize the production of natural cork stoppers, the most valuable product of the cork industry. Due to the very high within-plot variation in cork characteristics that is found on cork oak stands (Ferreira et al., 2000; Corona et al., 2005; Almeida and Tomé, 2010) it is difficult to make a decision on the optimum rotation period for cork extraction in a stand as the optimum is not the same for all the trees. This decision can only be achieved with the support of a growth and yield model that incorporates the prediction of individual tree cork caliber evolution over time. To date no such model is available in international literature. The development of a model to predict the evolution of cork caliber over time has to be based on the understanding of the process of cork formation which is well described in the cork oak related literature (e.g. Pereira and Tomé, 2004; Pereira, 2007; Paulo and Tomé, 2010). The first cork that the tree produces by the activity of the first periderm is called the virgin cork. Virgin cork is characterized by deep fractures and cracks that extend irregularly, but mostly longitudinally, due to the rapid radial growth of wood in young trees. The removal of the cork layer exposes the phellogen, causing this layer to consequently die. As a response, a new periderm (the traumatic periderm) is formed in order to protect the living tissues that were exposed to the environment. The unprotected tissues that dry out and die during this process fracture easily due to the radial growth of the new periderm that follows. These layers form the external part of the new cork layer, and are called the cork back. Internally to the cork back, seasonal cork rings develop until the next cork extraction takes place. The innermost cork layer is called the cork belly. The cork produced by this second periderm is called the second cork. It also displays deep fractures and is therefore still not the most desired for the production of cork stoppers. When the second cork is removed, the process of traumatic periderm formation is repeated and a new layer of cork is produced. The cork obtained from the third and subsequent cork extractions is named the mature or the reproduction cork. Cork oak is debarked during the growing season, starting the growth of new cork shortly after the extraction. This is the reason why the first and last growth rings are not complete. Fig. 1 (adapted from Natividade, 1940) shows an example of a cork with 9 years of age (8 complete growth rings). The first cork growth ring, located adjacent to the cork back, is usually smaller than the subsequent rings as it corresponds to the growth of cork in the growing season during which cork is extracted, after debarking. Therefore it is not a complete growth ring and is usually called the 1st half ring. Similarly, the last cork growth ring corresponds also to a partial growing period, before cork extraction, and is called the last

half ring. The production of new cork is more intense in the years immediately after debarking and decreases from then on (Pereira, 2007). The caliber of a cork plank with tc years (cork age is designated by tc in order to distinguish from stand/tree age that is usually designated by t) is equal to the sum of the thickness of the tc − 1 complete rings plus the thickness of the two half rings (initial and final) and of the cork back thickness. For sake of clarity, the term cork caliber (cc) is used here for the total cork thickness while the thickness of complete rings is designated by thickness of cork complete rings or simply by cork thickness (ct) as did other authors (Sánchez-González et al., 2008). There are just a few cork growth models available in the literature (Tomé et al., 1998, 1999, 2001; Almeida and Tomé, 2008; Sánchez-González et al., 2008), all of them referring to the complete rings, and based on limited data sets. The first published cork growth model was developed for cork oak stands in Portugal on the basis of 43 cork planks obtained from industry and therefore coming from different, non-identified sites. The model was based on a difference equation (multiple asymptotes) derived from the monomolecular growth function (Tomé et al., 1998, 1999). Another cork growth model was developed by Almeida and Tomé (2008) using data from 189 cork samples coming from different sites but not covering all the area of distribution of cork oak in Portugal. The difference equation method was again used in the modeling of cork growth and the selected model was a difference equation derived from the Lundqvist function with the asymptote as the free parameter. Almeida and Tomé (2008) focus the problem of the prediction of cork caliber from the corresponding cork thickness of complete rings, proposing a hyperbolic function for the prediction of the ratio cork caliber/cork thickness as a function of cork age. Sánchez-González et al. (2008) developed a cork growth model for a specific cork oak forest in Spain using the generalized algebraic difference approach (GADA), selecting a GADA formulation derived from the Richards model by Cieszewski (2004). The main objective of this work was the development of a general model for the prediction of the evolution of individual tree cork caliber over time, based on the measurements taken at a certain point in time (forest inventory). The model should be applicable to the whole country and appropriate to be included in a growth model for cork oak stands management, where it can be used to predict caliber, a regressor in the mature cork biomass prediction model of Paulo and Tomé (2010). When incorporated in such type of models, it will contribute to improve the economic analysis of the impact of decreasing or extending the interval between two consecutive cork extractions. It will also allow for the prediction of cork biomass growth, which is essential for the prediction of the evolution of carbon stocks. A second objective of this research was to define an index to express cork growth potential. Within stand cork growth variability is very large and cork growth potential is not strongly correlated with site characteristics (e.g. Sánchez-González et al., 2007). There is a need to find a way to express the cork growth potential of each tree so that the site potential may be defined through the distribution of the tree cork growth potential.

2. Material and methods 2.1. Developing the methodology to predict the evolution of individual tree mature cork caliber Forest inventory of cork oak stands focus on both tree and cork characteristics. These last are obtained from the analysis of 20 cm × 20 cm cork samples taken at breast height (commonly

A.M. Almeida et al. / Forest Ecology and Management 260 (2010) 1303–1314

1305

Table 1 Mean annual temperature (T) and mean annual precipitation (P) of the ‘concelhos’ that were sampled and indication of the respective number of sampling points. Id

‘Concelho’

No. of plots

T (◦ C)

P (mm)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

Macedo Cavaleiros Meda Idanha-a-Nova Abrantes Chamusca Ponte de Sôr Portalegre Avis Mora Coruche Benavente Palmela Montemor-o-Novo Évora Alcácer do Sal Portel Moura Grândola Santiago do Cacém Sines Odemira Silves São Brás de Alportel

1 1 1 1 1 2 2 1 1 3 1 1 3 2 2 1 1 2 2 1 2 1 1

14 11 15 16 16 15 16 16 16 16 16 16 16 16 16 17 18 16 16 17 16 17 17

731 637 663 715 737 680 811 637 650 657 553 615 690 649 587 644 484 668 663 695 701 680 820

referred as “calas” by the forest managers and landowners) in every tree inside the inventory plot or, most of the times, on a sub-sample of trees. Cork caliber of each cork sample is measured before and after boiling. The prediction of the evolution of cork caliber in a growth and yield model for cork oak stands, given the value of cork caliber measured for each tree at a certain point in time (during the forest inventory), implies the following stages, each one requiring the development of a model: 1. Model for the prediction of the thickness of cork complete rings (cork thickness) at initial age (tc0 ) given the values of cork caliber measured in the forest inventory and the respective age: cttc0 = f (cctc0 , other variables) where cttc0 and cctc0 are cork thickness and cork caliber measured at age tc0 . 2. Cork growth model (for complete rings) to update cork thickness of every tree till the end of the desired rotation period: cttc = f (cttc0 , tc0 − 1, tc − 1) where cttc and cttc0 are the thickness of complete cork rings at cork age tc and the cork thickness at an initial cork age tc0 , respectively. Note that this model applies just to the complete growth rings; reason why the model uses (tc0 − 1) and (tc − 1) as regressors. 3. Model for the prediction of cork caliber at age tc, given the corresponding value of the thickness of complete rings, predicted with the previous model: cctc = f −1 (cttc , other variables) where cctc and cttc are cork thickness and cork caliber estimated for cork age tc. 2.2. Data The data used to develop the cork growth model came from 34 sampling points (Table 1) distributed across different regions

Fig. 2. Distribution of Quercus suber L. in Portugal and location of the 23 ‘concelhos’ where the cork samples were taken. The numbers correspond to the ‘concelhos’ presented in Table 1.

representative of the area of expansion of cork oak in Portugal (Fig. 2). The location of the sampling points was based on the relative importance of each cork production region in terms of area occupied by cork oak. The sampling procedures are described by Pereira (1999). At each sampling point the 20 closest trees selected for debarking at the time of the forest inventory were sampled. Trees with sanitary or morphological problems (e.g. forked at a very low height and with many branches) or otherwise damaged were avoided. A cork piece of approximately 20 cm × 20 cm (“cala”) was taken from each sample tree at the diameter at breast height. Cork caliber before and after boiling was measured in each cork sample as well as the width of each one of the complete growth rings. The measurement of the cork rings took place after boiling the samples. Cork caliber after boiling is, on average, 12.6% greater than the corresponding pre-boiling caliber. This average value was obtained from all the samples of cork available in the database SUBER-DATA (Coelho et al., 2002). A total of 662 samples of mature cork, aged between 8 and 16 years, were used for cork growth modeling, corresponding to 5832 periods of growth. It was impossible to identify the cork growth rings in 18 samples that were discarded from the analysis. As mentioned in the description of the methodology, only complete growth rings were considered in the development of the cork growth model and the dependent variable was the thickness of complete cork rings, here designated simply by cork thickness (ct). Table 2 briefly characterizes the data used for fitting the cork growth models. The models to predict cork caliber and cork thickness (complete rings) were developed with the data set used in the cork growth modeling. Additional 40 cork samples were added to the

1306

A.M. Almeida et al. / Forest Ecology and Management 260 (2010) 1303–1314

Table 2 Descriptive statistics for several variables of the data set used to develop the cork growth model. Variable

n

Mean

sd

Min

Max

tc (years) cc (mm) ct1 (mm) ct2 ct3 ct4 ct5 ct6 ct7 ct8 ct9 ct10 ct11 ct12 ct13 ct14 ct15

662 662 662 662 662 662 662 662 662 653 321 103 51 41 18 6 5

10.0 33.9 4.0 8.0 11.7 14.9 17.7 20.4 23.1 25.6 27.4 29.9 33.6 34.4 39.6 34.0 33.9

1.4 8.6 1.4 2.6 3.7 4.5 5.2 5.9 6.5 7.0 7.5 7.8 9.1 8.7 9.7 8.4 8.0

8 16.5 0.9 2.0 3.3 4.6 5.5 6.9 9.7 11.2 12.2 15.8 17.5 18.7 22.1 23.7 25.5

16 67.2 8.7 18.3 25.4 28.7 32.7 37.4 41.3 45.3 50.8 54.0 58.8 52.7 55.7 44.4 44.3

tc – cork age; cc – cork caliber; cti – accumulated cork thickness in the first i complete annual rings; sd – standard deviation.

data set, which are representative of corks with age between 5 and 8 years, collected in the ‘concelho’ (small administrative division of the country) of Coruche. Table 3 summarizes the data set used in the prediction of cork caliber and cork thickness. 2.3. Model for cork growth 2.3.1. Modeling methodology and candidate models Cork growth data has a structure similar to stem analysis data which led previous authors (Tomé et al., 1998, 1999, 2001; Almeida and Tomé, 2008; Sánchez-González et al., 2008) to adopt site index curves development methodologies for cork growth modeling. The algebraic difference approach (ADA) as well as the generalized algebraic difference approach (GADA) were used in the present study. The ADA consists in the formulation of a growth function in such a way that one of the parameters is allowed to vary between sites (site-specific or free parameter) while the other parameters are assumed to be common to all sites. For cork growth models, the free parameter was allowed to vary not only between sites but also between trees within each site. The method is attributed to Bailey and Clutter (1974) and has been used by many other authors (e.g. Borders et al., 1984; Tomé, 1989; Amaro et al., 1998; Elfving and Kiviste, 1997; Carvalho and Parresol, 2005). When applied to cork growth, the resulting equations have the following general formulation: cttc = f (ct0 , tc0 − 1, tc − 1, ˇ), where cttc is the cork thickness (complete rings) at cork age tc, ct0 is the cork thickness at an initial cork age tc0 and ␤ is a vector of parameters. Note that the age regressors are (tc0 − 1) and (tc − 1) because the model refers to complete cork growth rings. The essence of GADA relies in the fact that it allows several parameters to be site-specific (Cieszewski and Bailey, 2000). In the case of cork growth, the within-plot variation is very high, implying Table 3 Descriptive statistics of several variables of the data set used to develop the model to predict cork caliber and cork thickness (complete growth rings). Variable

n

Mean

sd

Min

Max

tc (years) ct (mm) cc (mm)

702 702 702

9.7 27.3 33.4

1.4 7.9 8.5

5 9.2 14.0

16 58.8 67.2

tc – cork age; ct - cork thickness; cc – cork caliber; sd – standard deviation.

the need for site-tree-specific parameters. The method consists in: (i) choosing a base equation to model the variable of interest (ct in this case), (ii) deciding which parameters of the function are to be related with a theoretical measure of site and tree growth potential (X), (iii) solving the equation for X and then (iv) obtain the dynamic model with the form ct = f(tc − 1, tc0 − 1, ct0 , ˇ), by substituting the solution of X in the explicit equation ct = f(tc − 1, X, ˇ) for the initial conditions tc0 − 1 and ct0 . In this study we have considered both the situations of one parameter and two parameters of the base model related to X. When only one parameter of the base model is related to X, the GADA is equivalent to an ADA model. Table 4 lists the ADA and GADA models considered in modeling cork growth. The first group of dynamic equations was derived from the Lundqvist (1957) growth function. The second group was based on the Richards (1959) function and the last group on the log-logistic (log-transformed) function. The different formulations of the same function are indexed by the free parameters: the parameter a is an asymptote, b is a parameter related with the growth rate and c is a shape parameter. The seven ADA models derived from the Lundqvist (L a and L b) (Tomé, 1989), Richards (R a, R b and R c) (Tomé, 1989) and log-logistic (LL a and LL b) functions, are one-local-parameter models, and they are either anamorphic (L a, R a and LL a) by choosing the asymptote parameter a to be site-tree-specific or polymorphic with single asymptotes (L b, R b, R c and LL b) by choosing either parameter b or c as site-tree-specific. The model of McDill and Amateis (1992) corresponds to the model derived from the log-logistic function with the b parameter being site-tree-specific. GADA models L ab, L abb, R ac and R acc, polymorphic models with multiple asymptotes, were developed by Cieszewski (2004) from the Lundqvist and Richards base models. Models LL ab1 and LL ab2, also polymorphic models with multiple asymptotes, were derived by Cieszewski (2002) from the log-logistic. 2.3.2. Model fitting and test of regression assumptions The equations originated by ADA and GADA can be classified as self-referencing functions to give the indication that the dependent variable at a given age is estimated from the value of the same variable at another age. An effective technique for parameter estimation in self-referencing equations needs to account for the individual trends in the data (Cieszewski et al., 2000; Cieszewski, 2004). Cieszewski et al. (2000) propose the use of the dummy variable approach in which all the data points are used to produce residuals without the need of any arbitrary choice regarding measurement intervals. In this technique the global parameters are simultaneously estimated with the site-tree-specific parameters, resulting in the same parameter estimates regardless of the selected base age. The models were fitted in the following form: ctij = f (cti0 , tc0 − 1, tcij − 1, ˇ) where ctij is the cork thickness of tree i at cork age tcij , cti0 is the cork thickness of tree i at an initial age tc0 (base age, that was set as equal for all trees) and ␤ is a vector of global parameters. Note that the cork growth model applies just to the complete cork growth rings, reason why the variables (tcij − 1) and (tc0 − 1) appear as regressors instead of tc and tc0 . nIn the fitting procedure, the variable cti0 is substituted by ct d , a sum of terms containing a site-tree-specific paramei=1 i0 i ter cti0 and a dummy variable di for each tree within each site which is equal to 1 for the tree i and 0 otherwise (n is the total number of trees sampled). Any initial cork thickness corresponding to an arbitrarily selected initial age tc0 (tc0 = 0 not allowed) is used as starting value for cti0 . The sum of terms collapses into a single parameter unique for each tree that does not depend on the initial age selected as base age.

LL ab1 Cieszewski (2002)

LL ab2 Cieszewski (2002) a1 +X0 1+b1 X0 t −c

LL b McD McDill and Amateis (1992)

a1 +X0 1+b1 /X0 t −c

LL a

R acc Cieszewski (2004)

R ac Cieszewski (2004)

Y=

Y=

a 1−(1−a1 /Y0 )(t0 /t)c

Y=

Y = Y0 1+bt0−c

1− exp(−bt0 )

1+bt −c

 1− exp(−bt) (c1 +c2 /X0 )

Y = Y0 ) F0 = ln(1 − exp( − bt0 ))

1/2 4b1 Y0 t0−c ) )

X0 = Y0 (1 + bt0−c ) a=X

X0 =

X0 = 0.5(Y0 − a1 + ((Y0 − a1 ) +

b=X

a = a1 + X b = b1 /X

1−b1 Y0 t −c 0

Y0 −a1

X0 = a = a 1 + X b = b1 X

0

Y0 t −c

a1 −Y0

2

Log-logistic (logtransformed) Y = a/(1 + bt−c )

Y = exp(X0 ) (1 − exp(−bt)) X0 = 0.5((ln Y0 − c1 F0 ) + (((c1 F0 − ln Y0 ) ) − 4F0 ) a = exp(X) c = c1 + (1/X)

2

1/2

) F0 = ln(1 − exp( − bt0 )) X0 = c=X

)

2

X0 = b=X

1/c

(1− exp(−bt0 ))c

Richards (1959) Y = a(1 − exp( − bt))c

Y0

X0 = 0.5((ln Y0 − c1 F0 ) + (((ln Y0 − c1 F0 ) ) − 4c2 F0 )

X0 = a=X

+

+ + ln Y0 + F0 ) F0 =

+ ln Y0 + F0 ) F0 =

a = exp(X) c = c1 +c2 /X

R c Tomé (1989) ln(Y0 /a) ln(1− exp(−bt0 ))

X0 =

1/2 − ln(1−(Y0 /a) t0

X0 = a = exp(X) b = b1 + (1/X)

a = exp(X) b = b1 + (b2 /X)

2

((b1 t0−c ((b1 t0−c 0.5(b1 t0−c 0.5(b1 t0−c

0

X0 = b=X

t −c 0

− ln(Y0 /a)

exp(−bt −c )

Y0

Solution for X with initial conditions (Y0 , t0 )

X0 =

Tree specific parameters

a=X Lundqvist (1957) Y = a exp( − bt−c )

Base model

Table 4 Candidate models for the cork growth model.

(c1 +(1/X0 ))

R b Tomé (1989) ) )

Y = a(1 − (1 − (Y0 /a)

1−exp(−bt0 )

Y = Y0

Y = a(Y0 /a)(ln(1−exp(−bt)))/(ln(1−exp(−bt0 )))

R a Tomé (1989)

(1/c) t/t0 c

L abb Cieszewski (2004) Y = exp(X0 ) exp( − (b1 + (b2 /X0 ))t−c )

 1−exp(−bt) c

L ab Cieszewski (2004) Y = exp(X0 ) exp( − (b1 + (1/X0 ))t )

1/2 ln Y0 ) + 4t0−c ) 1/2 2 −c ln Y0 ) + 4b2 t0 )

Y = a(Y0 /a)

(t0 /t)c

−c

L b Tomé (1989)

Model Id Reference(s)

L a Tomé (1989)

Dynamic equations

Y = Y0 exp(b(t0−c − t −c ))

A.M. Almeida et al. / Forest Ecology and Management 260 (2010) 1303–1314

1307

All the models were fitted with the dummy variable approach implemented within the PROC MODEL procedure of the SAS statistical software (SAS Institute Inc., 2004). The regression assumptions of normality, homocedasticity and independence of errors were all evaluated for every tested model. The normality of errors was assessed with QQ plots and, when needed, the violation of the normality of model errors was overcome by applying the Huber function for robust estimation (Myers, 1986). Possible violations of the homocedasticity assumption, which would lead to the use of weighted regression, were examined by plotting standardized residuals versus predicted values of cork thickness. The presence of correlation among observations of the same tree (cork sample) was assessed through the plots of standardized residuals (rstd) versus the respective lagged standardized residuals (lagged residuals were considered as missing values if belonging to a different plot) and, if appropriate, the autocorrelation structure was modeled with an autoregressive model. Box-plots of residuals for the different plots were analyzed in order to detect large differences among plots that would indicate the need to take into account the hierarchical structure of the data. 2.3.3. Model evaluation and selection of the final model The fitting ability of the models tested was evaluated with the residual sum of squares (SSres) and the percentage of variance explained by the model, equivalent to the coefficient of determination for linear regression (R2 ). Validation of the models is a very important step in model evaluation because quality of fit does not necessarily reflect the quality of prediction. However the issue of which data set must be used for model validation as well as the methodology that must be used is controversial (e.g. Soares et al., 1995; Vanclay and Skovsgaard, 1997; Kozak and Kozak, 2003; Wang et al., 2008). A new independent data set is ideal for model validation, but most of the times there is no independent data set available as one wants to optimize the data set used for fitting. Since a validation data set was not available for the present research, validation was based on the PRESS residuals. A PRESS residual is equivalent to the residual that would be obtained for observation i with a model fitted with a data set from which this observation had been excluded. PRESS residuals are therefore independent from the data set used to fit the model and their use avoids data splitting. The designation PRESS comes from the PRESS (prediction sum of squares) statistic that is often computed with these residuals. Several statistics were calculated with the PRESS residuals to analyze the predictive ability of the models. The mean value of the PRESS residuals (Mrp) was used to evaluate the bias of the models. The mean of the respective absolute values (Marp) and the 95th and 5th percentiles of the PRESS residuals distribution were used to evaluate the precision of the models. The modeling efficiency, or proportion of variability explained by the model, was also computed using the PRESS residuals. The possible relationships of bias and precision with cork caliber or cork age were also evaluated by computing the Mrp and Marp statistics by cork caliber class or cork age (lack of fit statistics in the terminology of Kozak and Kozak, 2003), and plotting them as function of cork caliber class and cork age, respectively. The selection of the final cork growth model used all the information available from the previous analyses as well as the knowledge and empirical evidence on the biology of cork growth. 2.3.4. Defining individual tree cork growth potential Using the cork growth model it is possible to define an individual tree cork growth potential index as the predicted thickness of complete cork rings at a certain cork age (base age). The concept is similar to the site index concept but applicable on a tree basis instead of on a stand basis. Cork extraction occurs usually 9 years after the last extraction, which suggests the age of 9 years as the

1308

A.M. Almeida et al. / Forest Ecology and Management 260 (2010) 1303–1314

Fig. 4. Relationship between cork caliber (cc) and the corresponding cork thickness of complete rings (ct).

These two models, when used in conjunction with the cork growth model for complete years, allow for the prediction of the evolution of cork caliber over time. A hyperbolic model performs reasonably well when cork thickness is larger than a certain value but it does not give reasonable predictions for the first years of cork growth. Therefore this model form was not used in the present research. A preliminary graphical analysis between cork caliber and the corresponding cork thickness of complete rings revealed that there is a linear relationship between the two variables (Fig. 4). The following models were therefore used to predict cork caliber and cork thickness in complete years: Fig. 3. Above: relationship between the ratio cc/ct (cork caliber and the corresponding cork thickness of complete rings) and the cork thickness of complete rings (ct). Below: mean and standard deviation of the ratio cc/ct per cork thickness class.

“natural” base age in order to express individual tree cork growth potential. Note however that the cork growth model developed is base age invariant, therefore any other cork age can be used instead of 9 years. Individual tree cork growth index (cgi) is then defined as the cork thickness in complete years after 9 years of growth. The growth potential of a stand can be expressed by the distribution of cork growth index within the stand. 2.4. Model to predict cork caliber and cork thickness in complete years 2.4.1. Modeling methodology Almeida and Tomé (2008) developed a model to predict cork caliber and cork thickness in complete years based on the relationship that was found between the ratio cork caliber/cork thickness and cork thickness. This relationship was also present in the data set used in the present research (Fig. 3). The model proposed by these authors assumes that the ratio cork caliber/cork thickness decreases nonlinearly with cork thickness and uses an hyperbole to model this ratio using cork age (tc) as a surrogate for cork thickness: cc = ct

tc (equation to predict cork caliber) ˛ + ˇtc

ct = cc

˛ + ˇtc (equation to predict cork thickness) tc

cc = ˛ + ˇct ct =

cc − ˛ ˇ

where cc is the cork caliber, ct is the corresponding thickness of complete rings and ˛ and ˇ are parameters to be estimated. The expression of the parameters ˛ and/or ˇ as a linear function of cork age was tested. 2.4.2. Model fitting and test of regression assumptions The two models to predict cork caliber and cork thickness in complete years were simultaneously fitted. The system includes endogenous regressors, therefore two stages least squares (SAS/ETS MODEL procedure, SAS Institute Inc., 2004) were used in order to avoid simultaneous equations bias. Cork age and cork growth index were used as instruments. Regression assumptions were tested as explained for the cork growth model, except for the presence of correlation among observations, which was not present due to the structure of the data (one measurement per tree). 2.4.3. Model evaluation and selection of the final model The evaluation of the fitting and prediction ability of the models was based on the methodology used in the development of the cork growth model. In the evaluation of the possible relationship of bias and precision with cork caliber and cork age the smaller and higher classes of cork caliber and cork age were merged with the adjacent classes, as the number of observations was very small. In this way it was ensured that all classes of cork caliber and cork age had at least 147 and 49 observations, respectively.

A.M. Almeida et al. / Forest Ecology and Management 260 (2010) 1303–1314

1309

Fig. 5. Standardized residuals of model LL a (rstd) versus lagged rstd for one annual lag period (left column) and two annual lag periods (right column). Top: no autocorrelation structure; middle: first-order autoregressive error structure AR(1); low: second-order autoregressive error structure AR(2).

3. Results 3.1. Cork growth model When fitting the candidate models, the estimates of the parameter b1 of model L abb, parameter m of the model R b (c = 1/1 − m) and of the parameter a1 of model LL ab1 were found to be not significantly different from zero; therefore reduced forms of models L abb, R b and LL ab1 were used with b1 = 0, c = 1 and a1 = 0, respectively. It was impossible to meet the convergence criteria for models R c and LL ab2, even after increasing the number of maximum iterations, decreasing the convergence criteria (to reasonable limits), testing different fitting methods or changing the initial parameter values. Therefore these models were discarded from further analysis. Plots of the standardized residuals against the predicted values of cork thickness (data not shown) showed no evidence of heterocedasticity in all the models. The normal probability plots for the standardized residuals (data not shown), gave indication of non-normality of the errors for all the models tested. This lack of normality was overcome by applying the Huber function. The analysis of the independence of the residuals showed the presence of

correlation among observations of the same tree (cork sample) and the autocorrelation structure was modeled with an autoregressive model. The analysis of the plots of standardized residuals (rstd) versus lagged rstd suggested a second-order autocorrelation structure (AR(2)) that proved to be the appropriate autocorrelation structure for all of the models tested. Fig. 5 presents, as an example, the plots of rstd versus lagged rstd for model LL a fitted with no autocorrelation structure; first-order autoregressive error structure AR(1) and second-order autoregressive error structure AR(2) for 1 and 2 lagged years. The fit statistics, R2 and residual sum of squares (SSres) values of the 11 models are presented in Table 5. The value of the free parameter a (asymptote) is also presented. The range of asymptote estimates for the multiple asymptote models was assessed by computing this parameter values using cork thicknesses equal to 11.2 and 44.8 mm, both at 8 complete years (9 years old cork). Table 5 shows that all the models fit the data well, explaining about 99% of the total variation. The asymptote estimates for the models derived from the Lundqvist function were too high and not consistent with biological knowledge of the cork growth. On the other hand, the models derived from the Richards function originate asymptote values that seem too low. Comparing the ADA and GADA models

1310

A.M. Almeida et al. / Forest Ecology and Management 260 (2010) 1303–1314

Table 5 Estimates of the asymptote, R2 , residual sum of squares (SSres) and prediction ability statistics values for the candidate models for cork growth modeling. Statistics of prediction ability 2

Model Id

Asymptote value (a)

R

SSres

Mrp

Marp

P95

P5

SSrp

Model eff.

La Lb L ab L abb b1 = 0 Ra R b c=1 R ac R acc LL a LL b McD LL ab1 a1 = 0

445.18–1780.74 1287.84 447.55–1735.34 573.85–1409.34 23.18–92.73 72.90 24.11–90.63 24.12–90.54 33.93–135.74 108.73 50.40–114.07

0.995 0.995 0.995 0.995 0.995 0.995 0.995 0.995 0.995 0.995 0.995

2596.6 2640.2 2596.5 2599.5 2660.4 2761.7 2658.9 2659.0 2639.5 2721.5 2646.8

0.019 0.051 0.032 0.035 0.014 0.029 0.023 0.023 0.010 0.017 0.022

0.576 0.595 0.580 0.580 0.585 0.600 0.586 0.587 0.583 0.599 0.588

1.328 1.405 1.355 1.349 1.365 1.403 1.366 1.367 1.333 1.373 1.359

−1.175 −1.152 −1.175 −1.161 −1.168 −1.181 −1.166 −1.169 −1.182 −1.194 −1.178

3548.34 3853.16 3600.87 3589.09 3620.87 3892.18 3657.91 3665.25 3613.95 3878.88 3681.35

0.993 0.993 0.993 0.993 0.993 0.993 0.993 0.993 0.993 0.993 0.993

Notes: The range of the estimates for asymptote (a) in multiple asymptote models was assessed by computing the parameter values using cork thickness equal to 11.2 and 44.8 mm, both at 8 complete years. Mrp – mean of the PRESS residuals; Marp – mean of the absolute value of the PRESS residuals; P95 and P5 – 95th and 5th percentiles, respectively, of the distribution of the PRESS residuals; SSrp – sum of squares of the PRESS residuals; model eff. – model efficiency.

no large differences in performance could be found. The residual sum of squares (SSres) for the GADA models are, in general, lower than for the ADA models, being lower in the models derived from the Lundqvist function. Table 5 also presents the values of the prediction ability statistics for each one of the candidate models. The model efficiency is similar for all the models, while the sum of squares of the PRESS residuals is slightly lower for the GADA models. In terms of bias there are no large differences. However the models L b and L abb are the most biased and the LL a model is the less biased. Generally, the GADA models are more precise, but the differences are very small, and some of the ADA models such as the L a and the LL a are comparable to the best GADA models in terms of precision of the predictions.

As a conclusion, models R a, R ac, R acc, LL a, LL bMcD and LL ab1 are the most accurate, presenting similar bias and precision. Fig. 6 shows, for these models, the mean value of the PRESS residuals and the mean of the absolute values of the PRESS residuals as a function of cork caliber and cork age, showing that this relationship is similar for all the models. Predictive bias is higher for older corks. However, the range of bias values is generally very small, varying from −0.022 to 0.068 mm when considering the cork caliber class and −0.145 to 0.395 mm for cork age. In terms of precision, the behaviour of the six models is also very similar. Precision is lower for the thicker corks (>54 mm) and for corks with 13 years (Fig. 6). Based on the various criteria considered in the evaluation of the models, fit and prediction statistics and value of the asymptote, the

Fig. 6. The mean value of the PRESS residuals (Mrp) and the mean of the absolute values of the PRESS residuals (Marp) as function of cork caliber class and cork age.

A.M. Almeida et al. / Forest Ecology and Management 260 (2010) 1303–1314

1311

Table 6 R2 , residual sum of squares (SSres) and prediction ability statistics values for the model to predict cork caliber (cc) and cork thickness of complete growth rings (ct). Statistics of prediction ability 2

Model

R

SSres

Mrp

Marp

P95

P5

SSrp

Model eff.

cc = f(ct) ct = f(cc)

0.972 0.972

1386.1 1238.0

−2.178E−06 −6.657E−06

1.115 1.053

2.479 2.000

−2.116 −2.343

1390.25 1241.63

0.972 0.972

Mrp – mean of the PRESS residuals; Marp – mean of the absolute value of the PRESS residuals; P95 and P5 – 95th and 5th percentiles, respectively, of the distribution of the PRESS residuals; SSrp – sum of squares of the PRESS residuals; model eff. – model efficiency.

model LL a was selected to model cork growth: cttc = ct0

1 + 19.50(tc0 − 1)

−1.088

1 + 19.50(tc − 1)−1.088

(1)

where cttc is the cork thickness at cork age tc, ct0 is the cork thickness at an initial cork age tc0 . Note that this model applies only to cork complete rings, reason why the variables (tc − 1) and (tc0 − 1) appear as regressors instead of tc and tc0 . Model LL a is an anamorphic model with multiple asymptotes. The two parameters were estimated with small standard errors, 0.4970 and 0.0145, respectively for parameters b (bˆ = 19.50) and c (ˆc = −1.088), indicating small confidence intervals for parameter estimates and model predictions. 3.2. Model to predict cork caliber and cork thickness in complete rings When expressing the ˛ and ˇ parameters as linear functions of cork age, the associated parameters were not significantly different from zero. Therefore, the following models were obtained to predict cork caliber and cork thickness: cc = 4.572 + 1.058ct

(2)

cc − 4.572 ct = 1.058

(3)

where cc is the cork caliber and ct is the corresponding cork thickness of complete rings. Both models fit the data well and present reasonable values for the prediction ability statistics (Table 6). Plots of the standardized residuals against the predicted values (data not shown) showed no evidence of heterocedasticity in the two models. The normal probability plots for the standardized residuals (data not shown), showed that the residuals have distribution close to the normal. The observed and predicted means of the ratio cc/ct per cork thickness class almost match, indicating that the model fits the data well (Fig. 7). Fig. 8 shows the bias and precision obtained with the two models as a function of cork caliber and cork age. Generally the models are very precise and the bias is small. The bias of both models is symmetrical, higher for corks with thickness lower than 27 mm and older than 11 years, and is smaller for corks with ages equal or lower than 8 years. In terms of precision, the models behaviour is similar. The predictive precision slightly decreases with cork caliber and is lower for younger corks.

where cgi is the cork growth index (thickness of complete cork rings at 9 years of age), cttc is the cork thickness at tc years, cctc is the respective cork caliber. The prediction of cork growth potential for each tree within a stand or plot can be made from one measurement of cork caliber (cc) provided that cork age (tc) is known. It involves two steps: 1. Computing the complete cork rings thickness (ct) using Eq. (3); 2. Computing cork growth index (cgi) using Eq. (4). Fig. 9 shows the predicted cork growth curves for four classes (15, 25, 35 and 45) of cork growth index (cgi) and the observed cork growth curves for the trees inside a randomly selected plot. The amplitude of the cork growth index values observed in the plot is close to the range of values found in the complete data set (the same result holds for the other plots). The cork growth models developed by Almeida and Tomé (2008), Sánchez-González et al. (2008) are graphically compared with the model proposed in this research – model LL a. The curves of the models developed by Almeida and Tomé (2008) and SánchezGonzález are similar to the model developed in this research, although the data sets and the fitting methods used are different. 4. Discussion The result of this study is a set of models to predict the evolution of mature cork caliber of individual trees over time. That includes three sub-models: a model to predict complete cork rings thickness at the initial age (time of forest inventory) from cork caliber; a cork growth model, applicable to the cork complete rings, and a model to predict cork caliber at the end of the projection period from the corresponding thickness of complete rings. The algebraic difference approach (ADA) as well as the generalized algebraic difference approach (GADA) were used in modeling

3.3. Cork growth index The cork growth index (cgi), defined as the cork thickness (complete cork rings) after 9 years of growth (8 complete years), can be obtained from the cork growth model: cgi = ct9 = cttc

=

1 + 19.50(9 − 1)−1.088 1 + 19.50(tc − 1)−1.088

cctc − 4.572 1 + 19.50(9 − 1)−1.088 1.058 1 + 19.50(tc − 1)−1.088

(4)

Fig. 7. Smoothed curves for mean observed (- - -) and predicted (. . .) cc/ct ratios for each cork thickness class versus the corresponding central value. Individual observed values are also represented.

1312

A.M. Almeida et al. / Forest Ecology and Management 260 (2010) 1303–1314

Fig. 8. The mean value of the PRESS residuals (Mrp) and the mean of the absolute values of the PRESS residuals (Marp) per cork caliber and cork age classes, for the models cc = f(ct) and ct = f(cc).

cork growth. Several dynamic equations with one or two site-treespecific parameters were analyzed, but two of the models did not converge, even after increasing the number of maximum iterations, decreasing the convergence criteria (to reasonable limits), testing different fitting methods or changing the initial parameter values. All the models were fitted with the dummy variable approach. The models were evaluated by their fitting ability, their prediction performance using PRESS residuals and their behaviour screened with the available knowledge of the cork growth process. The performance of the 11 models for which the convergence criteria were achieved was very similar and no relevant differences were found between the ADA and GADA formulations of the growth functions. Sánchez-González et al. (2008) also used ADA and GADA

formulations of growth functions to model cork growth, but using the nested iterative method of the stochastic regression approach to fit the models instead of the dummy variable approach. These authors faced more problems in their work with the convergence of the models than in the present research (the convergence criteria was not met for 5 out of the 10 candidate functions that they tested). These authors also obtained very similar performances for the ADA and GADA formulations and selected a GADA formulation of the Richards function based on a slightly better value of the AIC value. The choice between a common or a multiple asymptote model is very important for site index curves modeling as the stands approach the asymptotic value of dominant height when they are usually harvested. It is expected that dominant height

Fig. 9. Left: predicted cork growth curves for four classes (15, 25, 35, 45) of cork growth index (cgi) and observed cork growth curves for the trees in a randomly selected plot. Right: comparison between the cork growth curves developed by Almeida and Tomé (2008), Sánchez-González et al. (2008) and model LL a (proposed in this work).

A.M. Almeida et al. / Forest Ecology and Management 260 (2010) 1303–1314 Table 7 Set of fitted equations to predict the evolution of individual tree cork calibre over time (1–4) and equation for prediction of pre-boiling cork caliber (5). 1. Model to predict cork growth in complete rings: cttc = ct0 2. Model to predict cork growth index: cgi = ct9 = cttc

1+19.50(tc0 −1)−1.088 1+19.50(tc−1)−1.088

1+19.50(9−1)−1.088 1+19.50(tc−1)−1.088

3. Model to predict cork thickness of complete rings, after boiling: ct =

cc−4.572 1.058

4. Model to predict cork caliber, after boiling: cc = 4.572 + 1.058ct 5. Prediction of pre-boiling cork caliber: ccpb =

cc 1.126

cttc – cork thickness (complete rings) with tc years (mm); ct0 – cork thickness with tc0 (initial age) years (mm); cgi – cork growth index; cc – cork caliber (mm); ct – cork thickness (mm) and ccpb – pre-boiling cork caliber in the tree (mm).

growth depends on the quality of the site and therefore that a multiple asymptote model must fit better the data as this type of model is consistent with the expected biological behaviour. Cork is extracted from the tree when it attains the size needed for the cork stoppers industry, between 27 and 40 mm, far away from the maximum values of cork thickness that have been observed in some rare cork samples that were taken from trees not debarked for a long time. These “old” cork samples can be very thick, therefore far away from the range of data available to fit the cork growth models. Sánchez-González (2006) report values of 110 mm (cork sample with 38 years) for Spain. In Portugal the largest value measured by the authors of the present research is equal to 180 mm. Maybe this is the reason why there is not a large difference in the performance of ADA and GADA models when fitted to cork growth data. Generally, the GADA models were slightly better but there were some ADA models in the set of models with better performance. Among these was the selected model LL a that was considered as the one that performed the best. It is an anamorphic model with multiple asymptotes. The two parameters were estimated with small standard errors of prediction, indicating small confidence intervals for the parameter estimates as well as for the model predictions. The models developed to predict cork caliber and cork thickness of complete rings were based in the linear relationship between the two variables and were fitted simultaneously through the two stages least squares procedure in order to avoid simultaneous equations bias. These models performed better, according to the fitting and prediction evaluation criteria that were used, than the ones proposed by Almeida and Tomé (2008). The models developed in this paper jointly with the model for prediction of pre-boiling cork caliber (Table 7) allow the simulation of the evolution of individual tree cork caliber during the cork rotation. This set of equations is therefore an important tool to assess the most favourable cork age for extraction in each stand in order to maximize the production of natural cork stoppers, the most valuable product of the cork industry. 5. Conclusions The objective of this work was the development of a model to predict the evolution of individual tree cork caliber over time, based on the measurements taken at a certain point in time, applicable to the whole country, and appropriate to be included as a module of a growth and yield model for cork oak stands. The set of equations developed in order to fulfill this objective are presented in Table 7 (Eqs (1)–(4)). These models, jointly with model 5 for prediction of pre-boiling cork caliber, will allow: • the prediction of complete cork rings thickness given the values of the corresponding cork caliber measured in the forest inventory; • the prediction of cork thickness (complete rings) with the growth model until the end of the desired rotation period;

1313

• the prediction of cork caliber given the predicted value of cork thickness (complete rings); • define the cork growth potential of each tree, using the cork growth index proposed; • to improve the analysis of the impact of decreasing or extending the interval between two consecutive cork extractions on cork production and respective economic input for the owner; • to predict one of the input variables of an existing model for prediction of mature cork biomass (Paulo and Tomé, 2010), which is essential for the prediction of the evolution of carbon stocks. Acknowledgements This work is part of the PhD research program of the first author which has been partially funded by the investigation supported by Fundac¸ão para a Ciência e a Tecnologia (FCT) through the scholarship (SFRH/BD/17733/2004). It has been supported by the projects CarbWoodCork (POCI/AGR/57279/2004) and MOTIVE (EU FP7-ENV-2008-1 226544). We acknowledge Helena Pereira for allowing the use of data collected under Project PAMAF 4053 “Caracterizac¸ão da qualidade tecnológica das cortic¸as portuguesas” and Augusta Costa for helping in cork growth data collection. References Adrados, J.R.G., Hernández, F.G., Haro, R.C., 2000. La predicción del calibre del corcho al final del turno y su aplicación al muestreo de la producción. INIA, Madrid. Almeida, A., Tomé, M., 2008. Sistema para predic¸ão do crescimento da cortic¸a. Silva Lusitana. 16 (1), 83–95. Almeida, A., Tomé, M., 2010. Field sampling of cork value before extraction in Portuguese ‘montados’. Agroforest. Syst. 79 (3), 419–430, Doi:10.1007/s10457009-9260-8. Amaro, A., Reed, D., Tomé, M., Themido, I., 1998. Modeling dominant height growth: eucalyptus plantations in Portugal. For. Sci. 44 (1), 37–46. Bailey, R.L., Clutter, J.L., 1974. Base-age invariant polymorphic site curves. For. Sci. 20 (2), 155–159. Borders, B.E., Bailey, R.L., Ware, K.D., 1984. Slash pine site index from a polymorphic model by joining (splining) nonpolynomial segments with an algebraic difference method. For. Sci. 30 (2), 411–423. Carvalho, J.P., Parresol, B.R., 2005. A site model for Pyrenean oak (Quercus pyrenaica) stands using a dynamic algebraic difference equation. Can. J. Forest Res. 35, 93–99. Cieszewski, C.J., 2002. Comparing fixed-and variable-base-age polymorphic site equations having single versus multiple asymptotes. For. Sci. 48 (1), 7–23. Cieszewski, C.J., 2004. GADA derivation of dynamic site equations with polymorphism and variable asymptotes from Richards, Weibull and other exponential functions. In: Second International Conference on Forest Measurements and Qualitative Methods and Management. University of Georgia, Athens, USA, pp. 248–261. Cieszewski, C.J., Bailey, R.L., 2000. Generalized algebraic difference approach: theory based derivation of dynamic site equations with polymorphism and variable asymptotes. For. Sci. 46, 116–126. Cieszewski, C.J., Harrison, M., Martin, S.W., 2000. Practical methods for estimating non-biased parameters in self-referencing growth and yield models. In: Daniel, B. (Ed.), Warnell School of Forest Resources. University of Georgia, Athens, GA, PMRC-TR-2000-7. Coelho, M.B., Godinho, J.M., Ribeiro, F., Freire, J., 2002. Suber data 2. Manual para o utilizador. Publicac¸ões GIMREF RTG 8/2002. Universidade Técnica de Lisboa, Instituto Superior de Agronomia, Lisboa, Portugal. Corona, P., Dettori, S., Filigheddu, M.R., Maetzke, F., Scotti, R., 2005. Site quality evaluation by classification tree: an application to cork quality in Sardinia. Eur. J. Forest Res. 124, 37–46. Elfving, B., Kiviste, A., 1997. Construction of site index equations for Pinus sylvestris L. using permanent plot data in Sweden. Forest Ecol. Manage. 98, 125–134. Ferreira, A., Lopes, F., Pereira H, 2000. Caractérisation de la croissance et de la qualité du liége dans une région de production. Ann. For. Sci. 57, 187–193. Kozak, A., Kozak, R., 2003. Does cross validation provide additional information in the evaluation of regression models? Can. J. Forest Res. 33, 976–987. Lundqvist, B., 1957. On the height growth in cultivated stands of pine an spruce 635 in Northern Sweden. Medd. fran Statens Skogforsk 47, 64 pp. McDill, M.E., Amateis, R.L., 1992. Measuring forest site quality using the parameters of a dimensionally compatible height growth function. For. Sci. 38 (2), 409–429. Myers, R., 1986. Classical and Modern Regression with Applications. Duxbury 640 Press, Boston, MA, 359 pp. Natividade, J.V., 1940. A propósito da idade das cortic¸as. Boletim da Junta Nacional da Cortic¸a, no 17. Junta Nacional da Cortic¸a, Lisboa, Portugal. Pereira, H., 1999. Caracterizac¸ão da qualidade tecnológica das cortic¸as portuguesas. Relatório do Projecto Pamaf 4053. Instituto Superior de Agronomia, 185 pp.

1314

A.M. Almeida et al. / Forest Ecology and Management 260 (2010) 1303–1314

Paulo, J.A., Tomé, M., 2010. Predicting mature cork biomass with t years of growth based in one measurement taken at any other age. Forest Ecol. Manage., Doi:10.1016/J.foreco.2010.02.010. Pereira, H., 2007. Cork: Biology, Production and Uses. Elsevier, Amsterdam. Pereira, H., Tomé, M., 2004. Cork oak and cork. In: Burley, J., Evans, J., Yongquist, J.A. (Eds.), Encyclopedia of Forest Sciences. Elsevier, pp. 613–620. Richards, F.J., 1959. A flexible growth function for empirical use. J. Exp. Bot. 10, 290–300. Sánchez-González, M., 2006. Modelo de crecimento y producción para monte alcornocal. Tesis Doctoral, Universidade Politécnica de Madrid. ˜ Sánchez-González, M., Calama, R., Canellas, I., Montero, G., 2007. Variables influencing cork thickness in spanish cork oak forests: a modelling approach. Ann. For. Sci. 64 (3), 301–312. ˜ Sánchez-González, M., Canellas, I., Montero, G., 2008. Base-age invariant cork growth model for Spanish cork oak (Quercus suber L.) forests. Eur. J. Forest Res. 127 (3), 173–182. SAS Institute Inc., 2004. SAS/ETS® 9.1 User’s Guide. SAS Institute Inc, Cary, NC. Soares, P., Tomé, M., Skovsgaard, J.P., Vanclay, J., 1995. Evaluating a growth model for forest management using continuous forest inventory data. Forest Ecol. Manage. 71, 251–265.

Tomé M., 1989. Modelacão do crescimento da árvore individual em povoamentos de Eucalyptus globulus Labill, Região Centro de Portugal. Ph.D. Dissertation, Technical University of Lisbon. Tomé, M., Coelho, M.B., Lopes, F., Pereira, H., 1998. Modelo de produc¸ão para o montado de sobro em Portugal. In: Pereira, H. (Ed.), Sobreiro e Cortic¸a”. Centro de Estudos Florestais, Lisboa, pp. 22–46. Tomé, M., Coelho, M.B., Pereira, H., Lopes, F., 1999. A management oriented growth and yield model for cork oak stands in Portugal. In: Amaro, A., Tomé, M. (Eds.), Empirical and Process-Based Models for Forest Tree and Stand Growth Simulation. Edic¸ões Salamandra, Novas Tecnologias, Lisboa, Portugal, pp. 271– 289. Tomé, M., Coelho, M.B., Almeida, A., Lopes, F., 2001. O modelo SUBER. Estrutura e equac¸ões utilizadas. Relatórios técnico-científicos do GIMREF RTG 2/2001. Centro de Estudos Florestais, Instituto Superior de Agronomia, Lisboa, Portugal. Vanclay, J.K., Skovsgaard, J.P., 1997. Evaluating forest growth models. Ecol. Model. 98 (1), 1–12. Wang, M., Borders, B.E., Zhao, D., 2008. Parameter estimation of base-age invariant site index models: which data structure to use? – Reply. For. Sci. 54 (2), 129–133.