Development and applications of a growth model for Pinus radiata D. Don plantations in El Bierzo (Spain)

A dynamic growth model for Pinus radiata D. Don plantations in El Bierzo (Spain) was developed with data from two inventories of permanent plots, of between 7 and 36 years old, established by the University of León. In this model, stand conditions at any point in time are defined by three state variables (stand basal area, number of trees per hectare and dominant height). The model includes three transition functions derived by the generalized algebraic difference approach to enable projection of the state variables at any particular time. Once they are known, the number of trees in each diameter class is estimated with a distribution function, by recovery of the parameters of the Weibull function by use of the moments method. Finally, a generalized height-diameter function and a taper function allow estimation of total or merchantable stand volume. The model provides satisfactory predictions for a time interval of three years. Simulation of the growth of four stands under two silvicultural regimes and two different sites confirm that the estimates provided by the overall model adequately represent the effects of both stand density and site quality. Other applications for the model are analysed and discussed.


Radiata pine in the area of study
Radiata pine (Pinus radiata D. Don) is an exotic species introduced in plantations worldwide because of its high growth rate and all-around advantages (e.g., Sutton, 1999).It is a well represented species in northern Spain, particularly in the Basque Country (160 000 ha) and Galicia (92 000 ha).
Although radiata pine was introduced in the region of El Bierzo relatively recently, it currently occupies an area of approximately 15 000 ha (Fernández-Manso et al., 2001).The oldest stands in the region were planted in the 1970s, with seeds initially obtained from nearby areas, such as Galicia.These stands were originally destined for production of mining timber, and were therefore planted at high densities and managed in short rotations (approximately 15 years); silvicultural practices were rarely carried out.Mining activities have now declined and trends are changing towards extended rotations and enhancement of silviculture practices.These afforestations introduced the species for the first time in a typical inland region of Spain, where Mediterranean and Eurosiberian conditions meet.
Most of the region is a rural marginal area in which forestry is becoming a highly recommended type of land use.Radiata plantations would certainly help to meet the challenges threatening the area (land abandonment, mine restructuring, decline of traditional economy and of population), thus contributing to sustainable development.

Forest growth models
Forest management decisions are based on information about current and likely future forest conditions.Consequently, it is often necessary to predict the changes in the system with growth and yield models, which estimate forest dynamics over time.Such models have been widely used in forest management because they enable updating of inventories, prediction of future yields, and exploration of management alternatives, thus providing information for decision-making in sustainable forest management (Vanclay, 1994).
The wide range of forest growth models available differ in complexity and the detail in which they describe the systems under consideration.At one end of the range are the traditional empirical models based on periodic tree measurements, which make no attempt to measure all factors that may affect tree growth, and at the other end the complex process models, based on the mechanisms inherent to growth, which incorporate a large amount of information in the response functions.
Of the two basic types of empirical models (wholestand and individual-tree models), whole-stand models are generally recommended when dealing with homogeneous, even-aged, pure stands (García, 1988(García, , 1993; Vanclay, 1994), because they can be constructed on the basis of variables often available in forest inventory data, and also represent a good compromise between generality and accuracy of the estimates.
Whole-stand models characterize the state of the stand by means of a small number of aggregate variables, such as basal area, mean diameter, volume per hectare, stems per hectare, average spacing or top height (García, 1993).This type of models therefore require few details for growth simulation.On the other hand, they provide rather limited information about the future stand (in some cases only stand volume) (Vanclay, 1994).To overcome underlying limitations, whole-stand models can be disaggregated mathematically using a diameter distribution function, which may be combined with a generalized height-diameter equation and with a taper function to estimate commercial volumes.Similar methodologies have been used by Mabvurira et al. (2002), Kotze (2003), Trincado et al. (2003), Diéguez-Aranda et al. (2006a) and Castedo-Dorado et al. (2007) in the development of forest growth models for plantations.
Preliminary studies in the region of El Bierzo have provided basic guidelines for the establishment and management of plantations (Castedo-Dorado et al., 2005a).However, the integrity and performance of these guidelines require a more elaborated simulation tool, such a dynamic growth model.
On the other hand, a recently developed whole-stand model for radiata pine stands is in use in the neighbouring region of Galicia (Castedo-Dorado et al., 2007).This model was developed with data from 225 inventories of permanent plots and it was demonstrated to be robust for medium term projections of stand volume.Nevertheless, it produces biased estimates when applied to the local data from the El Bierzo region, as explained in subsequent sections in this study.For this reason, an attempt was made to localize the already fitted model, so that it can be applied to El Bierzo without the need to develop a new model.Most adaptation procedures fall into one of three categories (Huang, 2002): the parameter re-estimation method, the proportional adjustment method or the regression adjustment method.According to Huang (2002), the most desirable approach for localizing any model is to refit the individual components (submodels) with local data.Therefore, the structure of the model adopted in this study was the same as that described by Castedo-Dorado et al. (2007).The same submodels were considered (and in most cases the same base equations), and a new set of parameters were reestimated from the local data.
The objective of the present study was to develop a management-oriented whole-stand model for simulating the growth of radiata pine plantations in El Bierzo.
The model is constituted by the same interrelated modules as those developed for the Galician model: a site quality system (both for dominant height growth and site index prediction), an equation for reducing tree number, a stand basal area growth function, and a disaggregation system composed of a diameter distribution function, a generalized height-diameter function and a total merchantable volume equation.

Data
The dataset used was obtained from two different sources.Initially, during the early winter of 2003, a network of 45 permanent sample plots was established by the University of León in pure radiata pine plantations in the area of study.The plots were located throughout the area of distribution of the species in the area of study, subjectively selected to cover adequately the existing range of ages, stand densities and sites.The plots were rectangular, and sized between 200 and 900 m 2 (mean size, 361.4 m 2 ), depending on stand density, to achieve a minimum of 50 trees per plot.Limits avoided stand border effects.A numbered label was nailed to the bark of every tree.
The measurements included the diameter at breast height (1.3 m above ground level) of each tree, which was measured twice (measurements made at right angles to each other) to the nearest 0.1 cm with callipers.The arithmetic mean of the two measurements was then calculated.Total height was measured to the nearest 0.1 m with a digital hypsometer (Vertex III) in a 30-trees randomized sample, and in an additional sample including the dominant trees.The number of domi-nant trees per plot was considered as the proportion of the 100 thickest trees per hectare (Assmann, 1970), thereby depending on plot size.Descriptive variables of each tree were also recorded (e.g. if they were alive or dead, or affected by disease).
A subset of 32 of the initially established plots was re-measured in the winter of 2006.These plots were selected for the dynamic components of the model.
The first source of data was the two inventories carried out in 2003 and 2006.The stand variables calculated for each plot and inventory were: age (t, years), stand basal area (G, m 2 ha -1 ), number of trees per hectare (N), square mean diameter (d g , cm), dominant height (H, m), defined as the mean height of the 100 thickest trees per hectare, and dominant diameter (D 0 ) defined as the mean diameter of the dominant trees.The stand age was estimated by counting rings from samples extracted with a Pressler borer from the base (stump-level) of three trees per plot and by direct knowledge of the age of plantation.
The mean, maximum and minimum values and the corresponding standard deviation for each of the main stand variables obtained from both inventories are shown in Table 1.
In addition, one or two dominant trees were destructively sampled at 23 locations.These trees were selected as the first two dominant trees found outside the plots but in the same plantations within ± 5 % of the mean diameter at 1.3 m above ground level and mean height of the dominant trees.The trees were felled to leave stumps of average height of 0.1 m; total bole length was measured to the nearest 0.1 m.The logs were cut at 1 or 2 m intervals.At each cross-sectional point, the diameter was measured and the number of rings was counted, and then converted to age above stump height.As cross-section lengths do not match the periodic

Development and fitting of transition functions for state variables
Transition functions must possess some important properties to provide consistent estimates: (i) consistency, i.e. no change for zero elapsed time; (ii) path-invariance, i.e. the result of projecting the state first from t 0 to t 1 , and then from t 1 to t 2 , must be the same as that of the one-step projection from t 0 to t 2 ; and (iii) causality, i.e. a change in status can only be influenced by inputs within the relevant time interval (García, 1994).Transition functions generated by integration of differential equations (or summation of difference equations when using discrete time) automatically satisfy these conditions.Fulfillment of the above mentioned properties for the transition functions depends on both the construction method and the mathematical function used to develop the model.Most of these properties can be achieved with techniques of dynamic equation derivation, known in forestry as the Algebraic Difference Approach (ADA) (Bailey and Clutter, 1974) or its generalization (GADA) (Cieszewski and Bailey, 2000).
Dynamic equations have the general form (omitting the vector of model parameters) Y = f (t, t 0 , Y 0 ), where Y is the value of the function at age t, and Y 0 is the reference variable defined as the value of the function at age t 0 .The ADA essentially involves replacing a base-model site-specific parameter with its initial-condition solution.The GADA allows expansion of the base equations according to various theories about growth characteristics (e.g., asymptote, growth rate), thereby allowing more than one parameter to be site-specific and allowing the derivation of more flexible dynamic equations (see Cieszewski, 2001Cieszewski, , 2002Cieszewski, , 2003)).More details about ADA and GADA derivation can be viewed in Cieszewski and Bailey (2000) and Cieszewski (2002).The ADA or GADA can be applied in modelling the growth of any site dependent variable involving the use of unobserv-height growth, Carmean's method (1972) was used to adjust height-age data from stem analysis to account for this bias.Log volumes were calculated with Smalian's formula, and the top of the tree was considered as a cone.Tree volume above stump height was aggregated from the corresponding log volumes and the volume of the top of the tree.The second source of data corresponds to the 41 felled trees, and allowed development of dominant height growth curves and total and merchantable volume equations.
Summary statistics, including mean, maximum, minimum and standard deviation of each of the main tree variables are shown in Table 2.

Model structure
The proposed whole-stand growth model is based on the state-space approach (Vanclay, 1994), which makes use of state variables to characterize the system at an initial stage, and transition functions to project all or some of the state variables in the future.In addition, it should be possible to estimate other variables of interest from the current values of the state variables through the so-called output functions (García, 1994(García, , 2003)).
The state of a system at any given time may be roughly defined as the information needed to determine the behaviour of the system from that time on; i.e., given the current state, the future does not depend on the past.According to García (1994), and considering that we are dealing with single-species stands derived from plantations in which different management regimes have been carried out, three state variables (dominant height, num-  able variables substituted by the self-referencing concept (Northway, 1985) of model definition (Cieszewski, 2004), such as stand height, number of trees per unit area, stand basal area, stand volume, stand biomass or stand carbon sequestration.
In the present study, dynamic equations with one and two site-specific parameters were tested.For projecting the state variables dominant height and stand basal area, we focused our efforts on the six dynamics equations used by Barrio et al. (2006Barrio et al. ( , 2008) ) for modelling the growth of these variables for maritime pine and poplar plantations.These equations are derived from three well-known growth functions (Korf -cited in Lundqvist, 1957-, Hossfeld -Hossfeld, 1882-and Bertalanffy-Richards -Bertalanffy, 1949, 1957;Richards, 1959) and have been widely used in modelling stand height and stand basal area growth (e.g.McDill and Amateis, 1992;Amaro et al., 1997;Falcao, 1997;Tomé et al., 2001;Cieszewski, 2002;Diéguez-Aranda et al., 2005).
The formulations of the equations analysed (M1-M6) are shown in Table 3.As general notational convention, a 1 , a 2 ,…, a n were used to denote parameters in base models, while b 1 , b 2 ,…, b m were used for global parameters in subsequent GADA formulations.
In addition, a dynamic equation was developed for predicting the reduction in tree number due to densitydependent mortality, which is mainly caused by competition for light, water and soil nutrients within a stand.Although many functions have been used to model empirical mortality equations, only biologically-based functions derived from differential equations include the set of properties that are essential in a mortality model (Clutter et al., 1983): consistency, path invariance and asymptotic limit of stocking approaching zero as old ages are reached.
In the present study, an equation for estimating reduction in tree number was developed on the basis of a differential function in which the relative rate of change in the number of stems is proportional to an exponential function of age: (1) where N is the number of trees per hectare at age t and α, β and δ are the model parameters.This function was selected by Álvarez et al. (2004) and Castedo-Dorado et al. (2007) to develop an equation in difference form for estimating reduction in stem number.In the present study, a dynamic equation was developed on the basis of this differential equation by use of the ADA, considering α as the site-specific parameter.
The individual trends observed in dominant height, number of trees per hectare and stand basal area data from the plots can be modelled by considering that individuals' responses all follow a similar functional form with parameters that vary among individuals (local parameters) and parameters that are common for all individuals (global parameters).Following Cieszewski et al. (2000), we estimated the random site-specific effects simultaneously with the fixed effects using a base-age invariant (BAI) parameter estimation technique: the dummy variables method.This technique has been used in several other studies in fitting growth functions for these stand variables (e.g., Diéguez-Aranda et al., 2005;Barrio et al., 2006Barrio et al., , 2008;;Castedo et al., 2007).
In the general formulation of the dynamic equations, the error terms are assumed to be independent and identically distributed with zero mean.Nevertheless, because of the longitudinal nature of the data sets used for model fitting, correlation between the residuals within the same individual (plot or tree) may be expected.This problem may be especially important in the development of the dominant height dynamic model on the basis of data from stem analysis, because of the number of measurements corresponding to the same tree.Nevertheless, in the construction of the dynamic equations for reduction in tree number and for basal area growth, which implies the use of data from the first and second inventory of 32 plots, the maximum number of possible time correlations among residuals is practically inexistent, and therefore the problem of autocorrelated errors can be ignored in the fitting process.
To overcome the possible autocorrelation, we modelled the error terms using a continuous time autoregressive error structure (CAR(x)), which allows the model to be applied to irregularly spaced, unbalanced data (Gregoire et al., 1995;Zimmerman and Núñez-Antón, 2001).To evaluate the presence of autocorrelation and the order of the CAR(x) to be used, graphs representing residuals plotted against lag-residuals from previous observations within each tree or plot were examined visually.The dummy variables method and the CAR(x) error structure was programmed by use of the SAS/ETS  MODEL procedure (SAS Institute Inc., 2004a), which allows for dynamic updating of the residuals.

Disaggregation system
Diameter distribution Among the parametric density functions that have been used to describe the diameter distribution of a stand (e.g., Charlier, Normal, Beta, Gamma, Johnson S B , Weibull), the Weibull function is the most frequently used, because of its flexibility and simplicity (Cao, 2004;Merganic and Sterba, 2006;Palahí et al., 2006).
Expression of the Weibull density function is: (2) where x is the random variable, a the location parameter defining the origin of the function, b the scale parameter and c the shape parameter controlling the skewness.
Of the two basic methodologies used to obtain the Weibull parameters (parameter estimation and parameter recovery), we used the parameter recovery approach, because as stated by several authors (e.g., Cao et al., 1982;Borders and Patterson, 1990;Torres-Rojo et al., 2000;Parresol, 2003) it provides better results than parameter estimation, even in long-term projections.To recover the Weibull parameters we used the momentsbased parameter recovery method (Burk and Newberry, 1984) because it directly warrants that the sum of the disaggregated basal area obtained by the Weibull function equals the stand basal area provided by an explicit growth function of this variable, resulting in numeric compatibility (e.g., Frazier, 1981).
In the moments method, the parameters of the Weibull function are recovered from the first three order moments of the diameter distribution (i.e. the mean, variance and skewness coefficient, respectively).Alternatively, the location parameter (a) may be set to zero.The use of this condition restricts the parameters of the Weibull function to two, therefore making it both easier to model and providing similar results to the three-parameter Weibull, at least for even-aged, single-species stands (Maltamo et al., 1995;Mabvurira et al., 2002).To recover parameters b and c the following expressions were used: where d is the arithmetic mean diameter of the observed distribution, var is its variance and the Γ is the gamma function.
Once the mean and the variance of the diameter distribution are known, while taking into account that Eq. 3 only depends on parameter c, the latter can be obtained by iterative procedures.Subsequently, parameter b can be calculated directly from Eq. 4. As the disaggregation system is developed for inclusion in a wholestand growth model, only the arithmetic mean diameter requires to be modelled, and the variance is directly obtainable from the arithmetic and the quadratic mean diameter (d g ) by use of the relationship var = d 2 g -d 2 .Hence, the arithmetic mean diameter was modelled with Eq. 5, which ensures that predictions of d are lower than d g : Where X is a vector of explanatory variables (e.g.dominant height, number of trees per hectare, age) that define the state of the stand at a specific point in time and must be obtained from any of the functions of the stand growth model and β is a vector of parameters to be estimated.
A diagram of the disaggregation system including all the components proposed in the present study is reported by Diéguez-Aranda et al., (2006 a ).

Height estimation for diameter classes
Once the diameter distribution is known, the next step is the estimation of the height of the average tree in each diameter class.Since the height-diameter relationship varies from stand to stand due to heterogeneous site conditions and silviculture state, and is not constant over time even within the same stand (Assmann, 1970), we used a generalized h-d model which takes into account stand variables that introduce the dynamics of each stand into the model (e.g., Curtis, 1967;Gadow and Hui, 1999).
The generalized h-d model used in the present study was a modification proposed by Castedo-Dorado et al. (2005b, 2006) on the basis of the Schnute (1981) function, by forcing it to pass through point (0, 1.3), to prevent negative height estimates for small trees, and to predict the dominant height of the stand (H) when the diameter at breast height of the subject tree (d) equals the dominant diameter of the stand (D 0 ).
Parameter estimates were obtained by ordinary least squares, by application of the Gauss-Newton's iterative method of the SAS/STAT ® NLIN procedure (SAS Institute Inc., 2004b).

Total and merchantable volume estimation
Once the diameter and height of the average tree in each diameter class are estimated, the total tree volume can be calculated directly by use of a volume equation.Since volume prediction to any merchantable limit is usually required, a taper equation can be used.Integration of a taper equation from the ground to any height provides an estimate of the merchantable volume to that height (Kozak, 2004).
Ideally, a volume estimation system should be compatible, i.e., the volume computed by integration of the taper equation from the ground to the top of the tree should be equal to that calculated by a total volume equation (Demaerschalk, 1972;Clutter, 1980).The total volume equation is preferred when classification of the products by merchantable sizes is not required, thereby simplifying the calculations and making the method more suitable for practical purposes.
A compatible system was fitted with data corresponding to diameter at different heights and total stem volume from 41 destructively sampled trees.To correct the inherent autocorrelation of the hierarchical data used, the error term was expanded by use of an autoregressive continuous model.The presence of autocorrelation and the order of the CAR(x) to be used were examined as explained in Section 2.2.2.
The compatible volume system of Fang et al. (2000) was used in this study because it was found to be very suitable for describing the stem profile and predicting stem volume for different species in several stand structures and regions (Diéguez-Aranda et al., 2006b;Corral et al., 2007), including Pinus radiata (Castedo-Dorado et al., 2007).The fittings were carried out by use of the SAS/ETS ® MODEL procedure (SAS Institute Inc., 2004a), which allows for dynamic updating of the residuals.
Aggregation of total (v) or merchantable (v i ) tree volume times number of trees in each diameter class provides total (V) or merchantable stand volume (V i ), respectively.

Selection of the best equation in each module
Comparison of the estimates of the different models fitted in each module was based on numerical and graphical analyses.Two statistical criteria obtained from the residuals were examined: the coefficient of determination for nonlinear regression (pseudo-R 2 ), showing the proportion of the total variance of the dependent variable explained by the model (Ryan, 1997), and the root mean square error (RMSE), which analyses the accuracy of the estimates.
Other important step in evaluating the models was graphical analysis of the residuals and examination of the appearance of the fitted curves overlaid on the trajectories of the dependent variables for each plot.Visual or graphical inspection is an essential point in selecting the most appropriate model because curve profiles may differ drastically, even though fit statistics and residuals are similar (e.g., Huang et al., 2003).

Overall evaluation of the whole-stand model
Although the behaviour of individual sub-models within the model plays an important role in determining the overall outcome, the validity of each individual component does not guarantee the validity of the overall outcome, which is usually considered more important in practice.The overall model outcome must therefore also be evaluated.The only method that can be regarded as "true" validation involves the use of a new independent data set (Vanclay and Skovsgaard, 1997;Kozak and Kozak, 2003;Yang et al., 2004).However, as new independent data for model validation were not available, only an evaluation of the predictive ability of the overall whole-stand model was carried out.For that purpose, observed state variables from the first inventory of the 32 plots measured twice were used to estimate total stand volume at the age of the second inventory, including all the components of the whole-stand model.It must be taken into account that total stand volume is the critical output variable of the whole model, since its estimation involves all the functions included and is important in economic assessments.
In order to assess whether the variance of the predictions is within some tolerance limits, the critical error statistic (E crit. ) was used.The critical error is expressed as a percentage of the observed mean and is computed by re-arranging Freese's χ 2 n statistic (Reynolds, 1984;Robinson and Froese, 2004): (6) where n is the total number of observations in the data set, is the observed value, ŷ i the value predicted from the fitted model, y is the average of the observed values, τ is a standard normal deviate at the specified probability level (τ = 1.96 for α = 0.05), and is obtained for α = 0.05 and n degrees of freedom.If the specified allowable error expressed as a percentage of the observed mean is within the limit of the critical error, the x 2 n test indicates that the model does not give satisfactory predictions; the contrary result indicates that the predictions are acceptable.
In addition, plots of observed against predicted values of stand volume were inspected.If a model is good, the slope of the regression line between observed and predicted values should pass through the origin at 45º.

Transition function for state variables
Among all the equations tested, model M4, the GADA form derived from the Hossfeld base equation was selected for height growth prediction and site classification: (7) where H 0 and t 0 represent the predictor dominant height (metres) and age (years) respectively, and H is the predicted dominant height at age t.
This function explained 99.35 % of the total variance of the data, and the value of the RMSE was 0.604 m.This model also showed adequate graphical behaviour, as it produced the most adequate height growth curves (Fig. 1, left-hand side) with a random pattern of residuals around zero with homogeneous variance and no discernable trend (Fig. 1, right-hand side).
As regards the transition function for reduction in tree number, a dynamic equation considering only one parameter to be site specific in the base model (Eq.8) described the data appropriately, providing the best results: (8) where N 0 and t 0 represent the predictor number of trees per hectare and age (years), and N is the predicted number of trees per hectare at age t.
Equation 8 explained 99.12 % of the total variance of the data and the RMSE was 57.02 trees/ha.The trajectories of observed and predicted number of trees over time for different initial density conditions are shown in Figure 2.
Regarding the transition function for stand basal area growth, among all the equations tested, model M5, derived from the Bertalanfy-Richards function solved by parameter a 2 provided the best results for the statistics used for comparison (R 2 = 0.982; RMSE = 1.74 m 2 ha -1 ), and provided a random pattern of residuals around zero (Figure 3).Moreover, the graphical analysis of the stand basal area growth curves overlaid on the trajectories of the observed values over time showed an adequate and biological behaviour.
The model is expressed as follows: ( ) where G 0 and t 0 represent the predictor stand basal area (m 2 ha -1 ) and age (years), and G is the predicted stand basal area at age t.
Application of the function for projecting stand basal area (G) requires an initial value of G at a given age, which may be obtained from diameter at breast height data.If this is not available, it must be estimated from other stand variables by use of an initialization equation.Following the methods of several authors (e.g., Amateis et al., 1995), we analysed several linear and nonlinear models with different explanatory stand variables (age, dominant height, site index, number of trees per hectare, relative spacing index, and combinations of these variables).Only data from 25 inventories, corresponding to ages younger than 15 years, were used, assuming that if projections based on ages older than this threshold are required, the initial stand basal area should be obtained directly from inventory data.A linear model with relative spacing index as independent variable behaved best.This model explained 76.4% of the total variance of the data, with a RMSE of 4.70 m 2 ha -1 .
The model obtained has the following form: G 0 = 53.97-1.318RS (10) where G 0 is the predicted stand basal area (m 2 ha -1 ) at a certain time, and RS is the relative spacing index in percentage.
The parameters of the stand basal area initialization function and of the three transition functions for the state variables were significant at a probability level of 0.05.

Disaggregation system
The equation selected for predicting arithmetic mean diameter and for use in the parameter recovery approach was: (11) where d is the predicted arithmetic mean diameter (cm), d g the quadratic mean diameter (cm) and t the stand age (years).The goodness of fit statistics were R 2 = 0.998 and RMSE = 0.168 cm.
The Kolmogorov-Smirnov test was used to compare the actual diameter distribution at the end of the projection interval with the distribution estimated by the moments method.The null-hypothesis tested was that the estimated distribution corresponds to the real distribution.In accordance with the results of similar studies, a significance level of 20% and a diameter class of one cm were considered.The results of the test showed that only one of the 32 diameter distributions estimated by the moments-based parameter recovery method was rejected.
Once the diameter distribution is known, the following generalized h-d relationship derived from the Schnute function was selected to estimate the height of the average tree in each diameter class: (12) where h is the predicted total height (m) of the subject tree, d its diameter at breast height (cm), and D 0 (cm) and H (m) are the mean diameter and mean height of the 100 thickest trees per hectare, respectively, of the stand where the subject tree is included.
This model showed a high predictive ability (R 2 = 0.934; RMSE = 1.474 m) and is very parsimonious because it only depends on two stand variables.
Finally, for total and merchantable volume estimation of the average tree in each diameter class, the compatible system based on the model developed by Fang et al. (2000) was fitted.This system is constituted by the following components: Taper function: (13) where I 1 = 1 if p 1 ≤ q i ≤ p 2 ; 0 otherwise I 2 = 1 if p 2 < q i ≤ 1; 0 otherwise p 1 and p 2 are relative heights from ground level where the two inflection points assumed in the model occur Merchantable volume equation: Volume equation: where: d = diameter at breast height over bark (cm); d i = top diameter at height h i over bark (cm); h = total tree height (m); h i = height above the ground to top diameter d i (m); h st = stump height (m); v = total tree volume over bark (m 3 ) above stump level; v i = merchantable volume over bark (m 3 ), the volume from stump level to a specified top diameter regression coefficients to be estimated; k = (π/40 000), metric constant to convert from diameter squared in cm in cm 2 to cross-section area in m 2 ; q i = h 1 / h.A third-order continuous autoregressive error structure was necessary to correct the inherent serial autocorrelation of the experimental stem data.The model provided a very good data fit, explaining 98.9% of the total variance of d i , and the RMSE = 0.847 cm.
All the parameters of the taper function and of all other equations corresponding to the disaggregation system were significant at a probability level of 0.05.

Overall evaluation of the whole-stand model
Assessment of model accuracy requirements was carried out by comparison between observed volume at the age of the second inventory and the volume estimated from the whole model.For this purpose, observed dominant height, number of trees per hectare and stand basal area from the first inventory of the 32 plots measured twice, served as initial values for the corresponding transition functions (Eqs.7, 8 and 9).These functions were used to project the stand state at the age of the second inventory.Equation 11 was then used to estimate the arithmetic mean diameter, which allowed calculation of the variance of the diameter distribution.Equations 3 and 4 were used to recover the Weibull parameters, thus allowing estimation of the number of trees in each diameter class.Finally, equations 12 and 15 were used to estimate the height and the total volume of the average tree in each diameter class respectively.Aggregation of total tree volume multiplied by the number of trees in each diameter class provided total stand volume.
A plot of observed against predicted values of stand volume obtained following the above procedure for the time interval considered (3 years) is shown in Figure 4 (left-hand side).The linear model fitted for the plot behaved well (R 2 = 0.983), however the plot showed that there was a slight tendency towards overestimation of stand volume for the prediction interval.A critical error of 10.79% was obtained in the projection of total stand volume for this time interval.
The plot of the observed stand volume values against the values predicted by the model developed with data from Galician stands (Castedo-Dorado et al., 2007) is shown in Figure 4 (right-hand side).Although the percentage variability explained by the linear model fitted for the scatter plot is high (R 2 = 0.978), the Galician model clearly underestimates the stand volume.Moreover, the critical error obtained in the projection of stand volume (23.84%) is much higher than that obtained by use of the local model.

Discussion
In this study, a dynamic whole-stand growth model for radiata pine plantations in an inland region of Spain (El Bierzo) is presented.The growth model described is comprehensive because it addresses all forest variables commonly incorporated in quantitative descriptions of forest growth.The method of construction adopted is robust because it is based on only three stand variables: dominant height (H), number of trees per hectare (N) and stand basal area (G).In accordance with García's state space approach to modelling plantations (García, 1983,1988), these variables summarize the historical events that affect stand development, and allow predictions from current state and future actions.The behaviour of the system is described by the rate of change of the state variables given by the corresponding transition functions.Besides, other stand variables (e.g.quadratic mean diameter, total or merchantable volume) can be obtained from the current values of the state variables.
In addition, the model can efficiently project stand development starting from different spacing conditions (e.g.densities) and considering different thinning schedules.
According to this structure, the whole-stand model only requires five stand-level inputs: age of the stand at the beginning and at the end of the projection interval (t 0 and t), and the initial dominant height (H 0 ), number of trees per unit of area (N 0 ) and stand basal area (G 0 ).Considering that modelling should be restricted to those variables available in forest inventory data (Burkhart, 2003), we consider our approach to be adequate for the practical management of P. radiata stands in El Bierzo, where predictor variables required are often measured and data are available at stand level.
The use of the GADA or ADA techniques in the development of the transition functions ensures that base-age and path invariance properties provide consistent predictions.Moreover, the functions were fitted by use of a base-age invariant (BAI) method that accounts for site-specific and global effects (Cieszewski et al., 2000).
The dominant height growth transition equation provided accurate values of site index from height and age and vice versa, regardless of the levels of site productivity.This transition function is one of the basic submodels in whole-stand growth models and allows assessment of the basic productivity on a stand-specific basis (e.g., Clutter et al., 1983).It is applied as follows: given site index (SI) and its associated base age (20 years in this case), to estimate the dominant height (H) of a stand for some desired age (t), simply substitute SI for H 0 and 20 for t 0 in equation 7. Similarly, to estimate site index at some chosen base age, given stand height and age, substitute SI for H and 20 years for t in the same equation.
Regarding the stand survival function, the methodology used for the projection estimates the reduction in the number of trees from a single mortality function fitted with all the data from plots measured in 2003 and in 2006.Although this interval of measurement collection is relatively short (Eid and Tuhus, 2001) it is considered representative of the natural mortality process that occurs in the stands, especially taking into account that radiata pine is a fast growing species.The accuracy of the survival function is important for generating realistic projections of the final output variables of the whole model (e.g., total or merchantable volume), especially when light thinnings or no thinnings are carried out (Avila and Burkhart, 1992), as was the case in most of the studied stands.
As regards the stand basal area projection equation, initial basal area and initial age provided sufficient information about the future trajectory of the basal area of the stand, regardless of initial spacing, thinning history or site quality.This can be explained because stand basal area and age are good enough estimators of quality in stands where silviculture is applied next to forest dynamics with low thinning (Bravo-Oviedo et al., 2004).The basal area initialization function developed should only be used to predict this variable in stands similar to those where the experimental data were collected, i.e., unthinned or lightly thinned stands younger than 15 years.The initialization and the projection function are not compatible because of the variation over time in the relative spacing index.This incompatibility is not a major problem, as the initialization function would only be used when no inventory data are available (Amateis et al., 1995).
Explanatory variables of the components of the disaggregation system can be easily obtained at any point in time from the three transition functions previously mentioned.The only exception is dominant diameter of the generalized h-d relationship, which is a difficult variable to project (Lappi, 1997), and must be estimated from the diameter distribution.
Total stand volume was selected in the present study as the critical output variable for the whole-stand growth model, although other stand variables can be assessed on the basis of this model.For instance, within the framework of climate change challenges, an increasing interest in studying the amount of biomass and carbon sequestration through forest management is expected and could be adequately tackled using the model.For example, if biomass equations that include the diameter at breast height and the total height as independent variables are developed for this species they can be easily incorporated into the disaggregation system proposed and converted in carbon pools (e.g., Merino et al., 2005).
Considering the required accuracy in forest growth modelling at stand-level, where a mean prediction error of the observed mean at 95% confidence intervals within ± 10%-20% is generally realistic and reasonable as a limit for the actual choice of acceptance and rejection levels (Huang et al., 2003), it can be stated that, on the basis of the critical error statistic obtained (10.79%), the model developed in this study provides satisfactory predictions.
Taking the same criterion into account, the model developed from the Galician data (Castedo-Dorado et al., 2007) cannot be applied directly to the region of El Bierzo, since the critical error statistics obtained in this case (23.84%) is higher than the upper limit usually considered for acceptance of a model.This result is expected a priori, since a fitted model is likely to produce errors when used on a local data set different from the data used in its construction (Huang, 2002).Moreover, these results also justify the need to adapt the Galician model, refitting the individual submodels on the basis of local (El Bierzo) data.
In summary, the overall evaluation of the model developed for El Bierzo demonstrates that it is robust, at least for short term projections of stand volume.This characteristic is especially interesting for planning purposes in fast growing species, as is the case here.As the study is based on stands of ages between 7-36 years, predictions for stands less than 7-years or more than 36 years-old should be used with caution.

Applications and conclusions
There are many potential applications for the model.As was explained, it can be applied with data normally available from stand inventories in the region.The model allows simulation of different types of thinnings (systematic, selective from below or semi-systematic).This may be extremely interesting, considering the need for implementation and generalization of silviculture practices in the area of study.The methodology used to simulate the thinning operations is that proposed by Alder (1979), which was successfully applied in subsequent studies (e.g., Álvarez et al., 2002).
As an example, Figure 5 shows the changes in basal area over the whole rotation for two silvicultural alternatives under two different site quality scenarios defined by site indices (SI) of 16 m and 26 m.The first alternative corresponds to the treatments usually applied in the past and occasionally at present: planting density (N 0 ) of 2500 stems per hectare (e.g., 2 × 2 m) without commercial thinnings, and rotation age of around 20 years.The second is a more intensive alternative, characterized by a wider initial space of 3 × 2.5 m (e.g., 1333 stems/ha), with two intermediate thinnings and an extended rotation age of 35 years.
The simulation given by the model seems realistic and sensitive both to initial stand density and site quality: dense plantations and better sites produce more stand basal area.
The model can also express yield in terms of product classes or log quality grade, which is vital management information for harvesting and marketing decision mak-ing (Vanclay, 1994).Information on size class distribution enables much more realistic evaluations of alternative silvicultural regimes in terms of both volume yield and financial returns.The predicted diameter distributions at clearfelling for the two silvicultural alternatives of the example are shown in Figure 6.As can be observed, the model takes into account the variation in the diameter distribution due to both stand density and site index: higher proportions of thicker trees are expected in alternatives with lower initial stand densities and on better sites.This is an example of how the use of the disaggregation system describes the stand much more thoroughly and readily allows simulation of thinning treatments.
The model also allows the development of curves of volume current annual increment (CAI) and mean annual increment (MAI).These curves are useful tools for the correct management of forest stands and contribute to estimation of the timing of intermediate and final cuts.As an example, both curves for the silvicultural schedule first considered, are shown in Figure 7.
It is well known that the interception between volume MAI and CAI curves indicates the biological rotation age for a given stand.The optimal rotation age when only average annual gain is considered is illustrated in Figure 7.As expected, biological rotation age is much shorter in better sites than in poor sites.
All the examples presented in this study can be directly used to optimise timber management planning and to evaluate alternative management regimes at stand level.These situations, in which planning can be carried out independently for each stand, will be the most typical in the region taking into account the typology of small non-industrial private forest owners.However, it is possible to run scenario simulations not only at the stand level but also at forest and regional level.In this case, where stable patterns of income are important, planning must be coordinated for all stands in the forest being considered (Clutter et al., 1983), and the information given by the model must be included in a forest level optimisation model.
In summary, the availability of the model can be considered important, taking into account that growth rates, together with management techniques, timber prices and financial returns are some of the most important factors driving forest management.In addition, the relatively simple structure of the growth model makes it suitable for embedding into landscape-level planning models and other decision support systems that enable  forest managers to generate optimal management strategies.
In conclusion, afforestation with radiata pine appears an adequate alternative considering both the stand growth rates that can be achieved and the socio-economic conditions of the area of study, in which most land would otherwise be abandoned.The high level of forest productivity confirms the results obtained by Romanyà and Vallejo (2004), who found a higher growth rate for this species in Mediterranean environments in Spain when compared to more Eurosiberian conditions.
The growth model developed can be considered as a potential tool for sustainable rural development.Multidisciplinary research should follow to complete radiata establishment, also considering wider objectives and derived environmental and socioeconomic synergies.Further consequences apply directly to rural development, socioeconomic and land management practices.

Figure 1 .Figure 2 .
Figure 1.Curves for site indices of 16, 20, 24 and 28 m at a reference age of 20 years overlaid on the trajectories of observed values over time (same tree or plot measurements joined by lines) (left-hand side), and pattern of residuals versus predicted values (right-hand side).

Figure 3 .
Figure 3. Stand basal area growth curves for stand basal areas of 15, 30, 45 and 60 m 2 ha -1 at 20 years overlaid on the trajectories of observed values over time.Same plot measurements joined by lines (left-hand side), and pattern of residuals versus predicted values (right-hand side).

Figure 4 .
Figure 4. Plots of observed stand volume values against model-predicted values, from the model developed in this study (lefthand side), and from the Galician model (Castedo-Dorado et al., 2007) (right-hand side).The solid line represents the linear model fitted to the scatter plot of data and the dashed line is the diagonal.R 2 is the coefficient of determination of the linear model.

Figure 5 .Figure 6 .
Figure 5. Changes in basal area over the rotation for the two silvicultural schedules simulated.

Figure 7 .
Figure 7. Volume mean annual increment (MAI) and current annual increment (CAI) for the silvicultural schedule corresponding to an initial stand density of 2500 stems per ha under two site conditions.

Table 1 .
Summarised data corresponding to the sample of plots used for model development

Table 2 .
Summarised data corresponding to the sample of 41 fallen trees used for model development

Table 3 .
Base models and ADA/GADA formulations considered for dominant height and stand basal area growth modelling