© INRA, EDP Sciences, 2010

1. INTRODUCTION

The European Beech (Fagus sylvatica L.) is one of the most important tree species in Europe. Its natural distribution range is defined by a maritime climate with moderate fluctuations of temperature and precipitation. The species avoids areas with extended very cold and dry periods and prefers moist sites. The natural range where beech is dominant and the phytogeographic range where it occurs, have been described in numerous research papers and more than 20 maps of the range have been published (Bolte et al., 2007).

In Swiss forests, European beech is the second species by growing stock with more than 18% (about 75 million m³) and the annual gross increments of this species between 1985 and 1995 represents the 20% of the total Swiss forests (SAEFL/WSL, 2005). This species will very possibly play an important role in Swiss forestry, because is the main natural species in lower altitudes where natural forests are politically demanded and it can produce highly valuable timber on good sites.

The yield tables developed by Paulsen (1795) possibly represent the first serious approach to beech growth modeling. These and subsequent generations of yield tables portray the development of pure beech stands in fixed (usually 5-yearly) time intervals for a given silvicultural treatment and growing site.

Due to improved computer technology and increasing information needs, the yield tables are being gradually replaced by flexible growth models. Several single tree growth simulators were developed during the past two decades to evaluate alternative silvicultural strategies (e.g. Nagel et al., 2002; Sterba and Monserud, 1997). Nevertheless, to date, the only growth model for the species in Switzerland is a yield table developed by Badoux (1983). This model provides limited information about the forest stand and does not reflect accurately the evolution under different stand density conditions. During the past years, there has been considerable interest in testing the applicability of individual tree models for a wide variety of Swiss forests in order to provide a potential alternative to the commonly used yield tables. However, for beech forests, the diameter increment was underestimated (Hallenbarter et al., 2005) and the lack of detailed stand structure information, not always available, had a strong effect on the quality of the simulation results (Schmid et al., 2006).

Our research has shown that it is possible to develop growth models which predict population mean parameters (e.g. dominant height) and area-based quantities (e.g. basal area per ha) with a comparatively high degree of accuracy. Thus, the objective of this contribution is to complement the range of existing beech growth models with a population-based approach, thus bridging the gap between yield tables and single tree models. We are using the extensive data base provided by the Swiss Federal Institute for Forest, Snow and Landscape Research to illustrate our state-space modeling approach.

2. MATERIAL AND METHODS

2.1. Data

The data used in this study originate form the growth and yield database of the Swiss Federal Institute for Forest, Snow and Landscape Research WSL in Birmensdorf Switzerland. The data were collected on 143 long-term pure and even-aged beech research plots (Fig. 1) on 103 locations, most of them on altitudes below 800 m in the Jura region, the Swiss plains and the Pre-Alps.

Figure 1 Map of all beech growth and yield plots (black dots = beech plots with tree height measurements, white dots = plots without height measurements). Map reproduced with permission of swisstopo (JA082265).

Some of the plots which are concentrated at the same location are thinning experiments where different thinning types were applied. In these experiments, the treatments included A-, B-, C-, D- and H-grade thinnings. These data reflect the response of beech stands of different ages growing on different sites to a variety of thinning regimes. The dataset used in this study covers a wide range of treatments, including extreme weights and types of removals. It is the best available source of information on beech growth in Switzerland.

The plots were re-measured between one and 17 times from 1889 to 2001. The shortest observation interval between plot remeasurements is two years, the longest 111 years. During the field enumeration, breast height diameters of all trees in the plot were measured crosswise with a caliper at a height of 1.3 m and the exact place of measurement was marked. Total heights and heights to crown base were assessed for selected trees including dominant trees, but not in all plots enumerations. Dominant height, H (m) was defined as arithmetic mean height of the 100 thickest trees per hectare and it was calculated for each plot where total heights were measured. In addition, descriptive variables such as the relative social position within the plot were assessed. Five of the plots are still being remeasured. The mean, maximum and minimum values and the standard deviation of the main stand variables are shown in Table I.

2.2. Modeling approach

The method used in this analysis is based on the state-space approach described by García (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 to future states.

2.2.1. State variables

In selecting the state variables, the principle of parsimony must be taken into account (García, 1994; Vanclay, 1995): the model should be the simplest one that describes the biological phenomena and that remains consistent with the structure and function of the actual biological system. A two dimensional vector including dominant height (H) and stand basal area (G) as explanatory variables may be sufficient to describe the state of the stand at a given age (t), and the advantage of including additional variables may not be justified (García, 1988). Nevertheless, in situations covering a wide range of silvicultural regimes, it may be necessary to include an additional variable such as the number of trees per unit area (N) (Castedo-Dorado et al., 2007a; Dieguez-Aranda et al., 2006; García, 1994). Taking into account that the effect of thinning on stand basal area growth was analyzed in this study, the stand conditions at any point in time are defined by three state variables (dominant height, number of trees per hectare and stand basal area) and the model include three transition functions. The outputs of these three transition functions provide the inputs for a static function, which predicts the total volume of the stand for a given age.

The three transition functions must possess some important properties to provide consistent estimates (García, 1994): (i) consistency, meaning no change for zero elapsed time; (ii) path-invariance, where the result of projecting the state first from t₀ to t₁, and then from t₁ to t₂, must be the same as that of the one-step projection from t₀ to t₂ and (iii) causality, in that a change in the state can only be influenced by inputs within the relevant time interval. Dynamic equations generated by integration of differential equations (or summation of difference equations when using discrete time) automatically satisfy these conditions.

Fulfilment of the above mentioned requirements depends on both the construction method and the specific mathematical function. A useful theory is the algebraic difference approach (ADA) developed by Bailey and Clutter (1974) and its generalization known as GADA (Cieszewski and Bailey, 2000). Dynamic equations have the general form (omitting the vector of model parameters) of Y = f(t, t₀, Y₀), where Y is the value of the function at age t, and Y₀ is the reference variable defined as the value of the function at age t₀. The ADA essentially involves replacing a base-model site-specific parameter with its initial-condition solution. The GADA allows expanding 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 derivation of more flexible dynamic equations (see Cieszewski, 2002; 2003).

In the present study, the transition functions for dominant height and stand basal area growth were derived with the generalized algebraic difference approach (GADA) and the transition function for reduction in tree number was derived with the algebraic difference approach (ADA).

An effective system for parameter estimation needs to be based on identifying individual trends represented in the data (Bailey and Clutter, 1974; Cieszewski, 2003). These trends represented in dominant height, number of trees per hectare and stand basal area can be modeled 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). In practice both base-age specific (BAS) and base-age invariant (BAI) estimation methods can be used. However, only BAI methods assume that the data measurements always contain measurement and environmental errors (both on the left and right hand sides of the model) that must be modeled (Cieszewski, 2003). Therefore, a BAI parameter estimation technique was used and from among the different BAI methods available the dummy variables method described in Cieszewski et al. (2000) was selected.

Based on the idea that one model parameter is specific for each individual (sample plot), and that the remaining parameters are shared by all individuals, Cieszewski (2003) proposed a simplified approach of mixed-effects model for transition functions development. The method of Cieszewski (2003), which we employed in this study, results from an adaptation of the nonlinear mixed-effects model definition by Lindstrom and Bates (1990), whose general formulation can be written as: Y_ij = f(A_iβ + B_ib_i, t_ij) + e_ij, where Y_ij is the jth response on the ith individual (i.e., dominant height or basal area of plot i at age measurement j); t_ij is the predictor for the jth response on the ith sample plot; f is a nonlinear function of the predictor and parameter vectors β (a p-vector of fixed population parameters –global parameters common to all sample plots) and b_i (a q-vector of the random effects associated with sample plot i – site-specific parameters); A_i is a design matrix of size r × p for the global parameters; B_i is a design matrix of size r × q for the site-specific parameters; r is the length of the vector of global parameters; and e_ij is a normally distributed random noise with zero mean.

Finding a model form that can describe the data well using only one parameter that varies with site is the first step of the method by Cieszewski (2003). If more than one parameter has to be site-specific, the varying parameters must be correlated with each other and defined as functions of fixed parameters and just one common varying parameter. After such a reparameterization, all except the varying parameter are uncorrelated parameters of the base model, and can be defined as pure fixed effects. The varying parameter is the unique parameter related to the random site effects (Lindstrom and Bates, 1990). The site-specific parameter of each sample plot can be defined as mixed-effects (both fixed and random effects) comprising the sum of the global mean for this parameter (ψ) plus a site-specific effect (ξ_i).

The next steps of the method by Cieszewski (2003) involve: (i) applying the GADA to ψ + ξ_i defined as X (a theoretical vector variable defining site quality with mean ψ and site-specific adjustments ξ_i with unknown distributional properties and an unknown magnitude of values); (ii) solving for X; (iii) substituting the solution for X in the explicit threedimensional site equation (Y = f(t,X)) for initial conditions t₀ and Y₀ (Y = f(t,t₀,Y₀), where Y₀ comprises for each sample plot the sum of the observed variable (dominant height or basal area) at age t₀ plus a site-specific effect e_0i interpreted as the random measurement and environmental errors); and (iv) fitting in one stage using a stochastic regression approach, or a varying parameter regression with desired fitting criteria. In the last fitting procedure, an estimation of the variables (i.e., Y₀ at t₀) resulting from fitting a curve to the observed trend in the data for each individual growth series are used as predictors instead of the observed variables. This direct use of data, as constants, does not violate regression assumptions because the measurement and environmental errors associated with these data are estimated simultaneously with all the other parameters of the model (Cieszewski, 2003).

As Cieszewski et al. (2000) pointed, the described base-age-invariant approaches are similar to mixed-effect modeling. The so-called local, or site-specific, parameters are in fact the random effects with a non-zero mean equal to the fixed effect representing average reference dominant height, basal area or number of stems per hectare at the given base age. This parameter is not modeled explicitly, but it can be easily recovered after the model fitting. All the other fixed effects are simply the global model parameters.

2.2.2. Transition functions for height and basal area

Dominant height and stand basal area growth models are key components of stand models, since they are directly related to economic variables such as stand volume. A number of researchers have pointed out desirable attributes of these two transition functions (Bailey and Clutter, 1974; Cieszewski and Bailey, 2000; Clutter et al., 1983). The most frequently listed criteria included the three general properties of transition functions described above and some specific requirements: polymorphism (for dominant height growth curves), sigmoid growth pattern with an inflexion point, horizontal asymptote at old ages (multiple asymptote for dominant height growth curves), logical behaviour (height should be zero at age zero and equal to site index at the reference age; the curve should never decrease), theoretical basis or interpretation of model parameters derived by analytically tractable algebraic operations and base-age invariance. The fulfilment of these attributes depends on both the model construction method and the statistical procedures used for model fitting, and cannot always be achieved.

Three different base equations were evaluated for modelling dominant height and stand basal area growth pattern: Bertalanffy (Bertalanffy, 1957); Log-logistic and Korf (1939). These functions have been widely used in modelling stand height and stand basal area growth (e.g., Álvarez González et al., 2005; Barrio Anta et al., 2006; Cieszewski, 2002; Diéguez-Aranda et al., 2006; Monserud, 1984).

As general notational convention, a₁, a₂,…, a_n were used to denote parameters in base models, while b₁, b₂,…, b_m were used for global parameters in subsequent GADA formulations. All the GADA based models have the general implicit form of Y = f (t, t₀, Y₀, b₁, b₂, ..., b_m).

The first GADA formulation was tested on the basis of the base model proposed by Bertalanffy (1957), whose integral form can be, for simplicity, represented as (1)where a₁ is an asymptote or limiting value, a₂ is often called a “rate parameter”, and a₃ is often referred to as “an initial pattern parameter”. To derive models with both polymorphism and variable asymptotes from Equation (1), more than one parameter has to be site-specific. Thus, both the asymptote a₁ and the shape parameter a₃ are assumed to be dependent on X. To facilitate such derivation the base equation is reparameterized into a form more suitable for manipulation of these two parameters (using exp (a′₁) instead of a₁ and taking the natural logarithm of the function) as follows: The site parameters can be conditioned to be consistently proportional to each other’s inverse over the site productivity dimension by the following definitions: which defines the following relationship: (2)For this formulation, the solution for X involves finding roots of a quadratic equation and selection of the most appropriate root to substitute into the dynamic equation. The selection of the most appropriate expression for X may depend on the equation parameters that in turn depend on the data and the domain of the applicable ages (Cieszewski and Bailey, 2000). The solution for X in Equation (2) with initial condition values t₀ and Y₀ is where Selecting the root more likely to be real and positive (i.e., the one involving addition rather than subtraction of the square root), and substituting it into Equation (2) results in the first GADA formulation used in this study: (E1)The next two GADA formulations were derived on the base of the log-logistic model: (3)where a₁ is the asymptote parameter, a₂ is the half-saturation parameter that defines the value of exp( − u) at which Y(u) = a₁ / 2, and u is usually a linear combination of logtransformations of the independent variable(s) (Monserud, 1984; Cieszewski, 2002).

Within the different forms investigated, the GADA based formulation by Cieszewski (2002) was found to perform particularly well. It considers the above-mentioned u as a logarithmically transformed function of age, such that u = a₃lnt. Therefore, Equation (3) can be simplified as: (4)In a first approach, Cieszewski (2002) replaced a₁ with a constant plus the unobserved site variable X, and a₂ by b₂ / X, becoming the next formulation (5)The solution for X involves again finding roots of a quadratic equation. With initial condition values t₀ and Y₀ it is Selecting the root more likely to be real and positive (i.e., the one involving addition rather than subtraction of the square root), and substituting it into Equation (5) results in the second GADA formulation: (E2)A second GADA formulation can be obtained based on Equation (4) by replacing a₁ with a constant plus the unobserved site variable X, and a₂ by b₂X(6)The solution for X in Equation (6) with initial condition values t₀ and Y₀ is Substituting it into Equation (6) results in the third GADA formulation: (E3)The fourth transition function is based on the base model proposed by Korf (1939): where a₁ is the asymptote parameter, and the other two parameters (a₂ and a₃) are related to the inflexion point and the growth rate.

Cieszewski (2002) replaced a₁ with the exponential of the unobserved site variable X, and a₂ by b₁ / X(7)The solution for X involves taking the natural logarithm of the functions and again finding roots of a quadratic equation. With initial condition values t₀ and Y₀ it is Selecting the root more likely to be real and positive (i.e., the one involving addition rather than subtraction of the square root), and substituting it into equation (7) results in the fourth GADA formulation: (E4)Different methods have been used to take into account the thinning effects in basal area growth. Some authors (Barrio et al., 2006; Knoebel et al., 1986; Zarnoch et al., 1991) have developed different equations for thinned and unthinned stands that have the same mathematical structure but that have been parameterized using different data sets. In the other hand, other studies (e.g. Hasenauer et al., 1997; Hynynen, 1995) have analyzed the inclusion of a thinning response function to express the basal area growth of a thinned stand as a product of a reference growth and the thinning response function: the reference growth accounts for the factors affecting stand growth in unthinned stands while the thinning response function predicts the relative growth response following thinning. In this work, we used the first approach to take into account the effect of thinning on stand basal area growth.

2.2.3. Transition function for reduction in number of trees

One of the most important transition functions of a dynamic growth and yield model is a mortality model that estimates the natural decline in number of trees caused by stand density and other environmental factors. However, mortality remains one of the least understood components of natural processes growth. Many functions have been used to model empirical mortality equations. Some of them are mathematical relationships between stem number reduction and stand variables (e.g. Álvarez González et al., 2004; Rennolls and Peace, 1986), others are biologically based functions derived from differential equations. These biologically based functions have properties that are essential in a mortality model but are not always present in a pure mathematical model (Clutter et al., 1983): consistency, path-invariance and asymptotic limit of stocking approaching zero as age becomes very large. Also, for even-aged stands it is reasonable to assume that ingrowth is negligible, so if age at time two is greater than age at time one, the density at time two will be less than density at time one. Specification of a difference equation model that possesses these properties and properly expresses the population relationship is often a difficult task (Clutter et al., 1983). In this work, the equation for estimating reduction in tree number was developed on the basis of the following differential function: (8)(9)(10)where N is the number of trees per hectare present at age t, ΔN / Δt is the instantaneous mortality rate operating at age t, and α, β, γ and δ the model parameters that regulate the mortality rate.

The three equations (Eqs. (8)–(10)) consider the mortality rate related to number of trees per hectare and age. Five different transition functions derived from these three differential functions have been used in this study, which are shown in Table II.

2.2.4. Model fitting

The state-space approach has been widely used to develop stand growth models (e.g. Castedo-Dorado et al., 2007a; Diéguez-Aranda et al., 2006; Nord-Larsen and Johannsen, 2007). Usually, each transition function is fitted separately, however the transition functions for dominant height, basal area and number of trees reduction together define a system of three regression equations and the random errors of the equations could be correlated. In order to improve the efficiency of the estimation the iterative seemingly unrelated regression (ISUR) method, that takes into account the cross-equation correlations, was used. The estimate of the cross-equation error covariance matrix to initiate the iterative procedure was obtained using Ordinary Least Squares. A similar approach was used by Nord-Larsen and Johannsen (2007) to develop a stand growth model for European beech in Denmark.

In the general formulation of the dynamic equations, the error terms e_ij 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 plot may be expected. To overcome any possible autocorrelation, we modelled the error terms using a 1-order continuous time autoregressive error structure (CAR(1)), which allows the model to be applied to irregularly spaced, unbalanced data (Gregoire et al., 1995): where e_ij is the jth ordinary residual on the ith sample plot, d₁ = 1 for j > 1 and d₁ = 0 for j = 1, e_{ij − 1} is the j − 1th ordinary residual on the ith sample plot, ρ₁ is the autoregressive parameter to be estimated, and t_ij − t_{ij − 1} is the time that separates the jth from the j − 1th observations within each sample plot. In such cases ε_ij now includes the error terms under conditions of independence.

To evaluate the presence of autocorrelation, graphs representing residuals plotted against lag-residuals from previous observations within each plot were examined visually. The heteroscedasticity was also examined visually, by plotting residuals as a function of predicted values, in order to investigate possible weighting factors.

The fitting of the equations were accomplished by SAS/ETS^® MODEL procedure (SAS Institute Inc., 2004). The correction for autocorrelation was programmed in the MODEL procedure by adding to the model formulation the term di*(rho**dist)*zlag1(resid.V), where di is a dummy variable which values is 0 for the first observation of each sample plot and 1 in other case; rho is the autoregressive parameter (ρ₁); dist is the time that separates the the jth from the j − 1th observations within each sample plot and zlag1(resid.V) is a SAS function referring to the previous residual of the stand variable V (dominant height, basal area or number of trees per hectare).

2.2.5. Model comparison and evaluation

The comparison of the estimates for the different equations fitted for each stand variable was based on numerical and graphical analyses of the residuals (e_ij). Three statistical criteria obtained from them were examined: bias (), which tests the systematic deviation of the model from the observations; root mean square error (RMSE), which analyses the accuracy of the estimates; and the model efficiency (MEF), which shows the proportion of the total variance that is explained by the model, adjusted for the number of model parameters and the number of observations. The expressions of these statistics are as follows: where y_i, and are the observed, predicted and average values of the dependent variable, respectively; n is the total number of observations used to fit the function; and p is the number of model parameters.

Other important steps in evaluating the functions fitted were 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 (Huang et al., 2003).

Once the complete model has been constructed, assessment of their validity is often needed to ensure that the predictions represent the most likely outcome in the real world (Yang et al., 2004). Moreover, although the behaviour of the functions for each stand variable 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 (Pretzsch et al., 2002; Yang et al., 2004) but the scarcity of such data forces the use of alternative approaches. The common method of splitting the data set in two portions does not provide additional information (Huang et al., 2003) and neither is recommended from the point of view of parameter estimation.

Moreover, other techniques such as double cross-validation or statistical testing provide very limited information about the predictive ability of the models (Kozak and Kozak, 2003; Yang et al., 2004). Therefore, we decided to defer model validation until a new data set is available for assessing the quality of the predictions.

The functions of each stand variable were evaluated to confirm their usefulness and reliability by the goodness-of-fit statistics and the graphical analysis described above (refer also to Fig. 2). In addition, plots of observed against predicted values of each state variable were inspected.

Figure 2 H-Residuals versus H-Lag1-Residuals for model E1 (first row), G-Residuals versus G-Lag1-Residuals for model E1 (second row) and N-Residuals versus N-Lag1-Residuals for model E5 (third row) fitted using a first-order continuous time autoregressive error structure.

To evaluate the maximum time interval for accurate projections, we used the critical error (E_crit). In this study, the critical error is expressed as a percentage of the observed mean basal area and is computed by re-arranging Freese’s statistic (Reynolds, 1984): where n is the total number of observations in the data set, y_i is the observed basal area, the basal area predicted from the fitted model, is the average of the observed basla area, τ 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.

3. RESULTS AND DISCUSSION

At first, the transition functions for H, G and N were fitted separately to select the best equation. The models were fitted using nonlinear least squares without expanding the error terms to account for autocorrelation. However, a slight trend in residuals as a function of lag1 residuals within the same sample plot was apparent in all the models. After applying a first-order continuous time autoregressive error structure, all the parameters associated with the error term were significant at the 5% level, and the error trends in residuals disappeared (Fig. 2). The sole purpose of the autocorrelation correction was to improve the interpretation of the statistical properties of the models, and it has no use in practical applications unless one is working repeatedly on the same sample plot.

The parameter estimates for each model, including its standard errors and p values, and the corresponding goodness-of-fit statistics are shown in Tables III to V for dominant height, basal area and reduction in number of trees, respectively. After considering the statistical criteria and analysing the graphical output, model E4, was selected as the best transition function for dominant height and basal area and model E5 for reduction in number of trees.

Although model E4, derived from the Korf function, does not achieve the best values based on the statistical criteria which were used to compare the transition functions for dominant height (Tab. III), in the graphical analysis this model showed the best performance.

For basal area growth, all the models explained more than 98% of the total variance, and model E4, derived from the Bertalanffy equation obtained the best values of the statistical criteria used (Tab. IV). However, as in the case of dominant height growth, the model E4, derived from the Korf equation showed a better performance.

Differences in basal area growth between thinned and unthinned sample plots were detected, therefore, different transition functions are proposed for each treatment. These results suggest that, for our data set, the stand basal area growth rate following a thinning exceeded the stand basal area growth rate of a stand with similar initial stand conditions which had not been recently thinned. This “thinning effect”agrees with several other studies (e.g., Amateis, 2000; Hasenauer et al., 1997; Pienaar and Shiver, 1984; Pienaar et al., 1985; Sharma et al., 2006). These studies involving different species have demonstrated that there is a difference, however slight, in the basal area growth in thinned and unthinned stands of the same initial age, site index and stand basal area.

This result, however, contradicts the findings of other studies. Barrio et al. (2006) and Castedo-Dorado et al. (2007b) proposed the same equation for thinned and unthinned stands assuming that the thinning effect was an inherent part of the model. The authors concluded that any contradictory results relating to the “thinning effects theory”may be at least partly attributed to the specific structure of the experimental data sets. Some data sets include only a narrow range of thinning trials where mainly low or moderate thinnings were carried out.

Due to the fact that the data used in this study do not provide enough information about the time interval during which the thinning effect is relevant for basal area growth, we propose that the transition function for basal area in thinned stands should not be used for projection intervals longer than 5 years.

Table V shows the parameter estimates for each of the 5 equations for the reduction in stem number and the statistics used to compare them. Models E5 to E7 are derived from a differential function considering that the relative rate of change in the stem number is proportional to a power function of age. Models E8 and E9 are derived from differential equations considering a hyperbolic and an exponential function of age, respectively.

Taking into account the statistics used to compare the models, there is no clear evidence of how to include age as an explanatory variable. However, only for model E5 all the parameters were significant. Therefore, this model was selected for fitting all the equations simultaneously. This equation implies that the relative rate of change in the number of trees is proportional to a power function of age. In the past, similar results were obtained for different species (e.g. Clutter and Jones, 1980; Dieguez-Aranda et al., 2005; Woollons, 1998).

Figure 3 Upper left: dominant height curves for site indices of 18, 22, 26, 30 and 34 m at a reference age of 80 years for the GADA equation E1. Upper right: reduction in tree number curves for 400, 300, 2 200, 3 100, 4 000 trees per ha at a reference age of 40 years for the model (E5). Down: basal area growth rates for thinned (left) and unthinned stands (right) as affected by initial stem number and age overlaid on the trajectories of basal area.

In a second approach, the selected models for H, G and N were fitted simultaneously to minimize the total sum of square errors in the system and to take the cross-equation error correlations into account. When ordinary least squares was used to fit the models simultaneously, the correlation value of the errors was only significant for the dominant height-basal area models (ρ = 0.3252). This value was reduced by 21.7% using iterative seemingly unrelated regression. The graphical analysis of residual H, G and N against the corresponding predicted values showed no trends.

The transition function for dominant height proposed for beech stands finally obtained by the simultaneous fitting is the following: where Y is the predicted dominant height (m) at age t (years), and Y₀ is the predictor dominant height at age t₀.

This model explained 98.1.5% of the total variance of the data, and its RMSE was 1.01 m. In selecting the optimum base age, it was found that 80 years was superior for predicting height for other ages. The curves for dominant height development for 18, 22, 26, 30, 34 m at a reference age of 80 years overlaid on the profile plots of the data set are shown in Figure 3.

The proposed transition functions for basal area for unthinned and thinned beech stands, respectively, obtained by the simultaneous fitting are the following: where Y is the predicted basal area (m²/ha) at age t (years), and Y₀ is the predictor basal area at age t₀.

The basal area transition function accounted for more than 98% of the total variance of the data with a RMSE of 1.27 m²/ha. The curves of basal area growth rates for thinned and unthinned stands as affected by initial basal area and age, superimposed over the trajectories of observed values over time are shown in Figure 3.

Finally, the proposed equation for estimating the natural stem number decline for beech stands is: where N is the predicted number of trees per hectare at age t (years), and N₀ is the predictor number of trees per hectare at age t₀.

The above model explained approximately 96.4% of the total variance of the data and the RMSE was 83.4 trees/ha. The trajectories of observed and predicted number of trees over time for 400, 1 300, 2 200, 3 100 and 4 000 trees per hectare at a reference age of 40 years are shown in Figure 3.

To achieve a satisfactory overall evaluation of the model, the observed basal area from any one inventory of the sample plots was projected forward for different time intervals using the transition function. All possible ascending growth intervals were considered and the thinnings were simulated at the same age that was carried out using the real rate G_removed / G_total. The root mean square error and critical error in basal area estimations were calculated for different time intervals and the results are presented in Figure 4.

Figure 4 Plots of RMSE in percent (up) and critical error in percent (down) for basal area predictions against different prediction intervals. The rectangle indicates the maximum time interval for a required accuracy of 15% in basal area estimations.

According to Huang et al. (2003) an accuracy within 10–20% of the observed mean value at 95% confidence intervals is reasonable as limit for the acceptance level. The 15% error is reached at a projection interval of 25 years and 18 years for RMSE and critical error, respectively. Therefore, we conclude that the model provides satisfactory predictions for time intervals up to 20 years.

4. CONCLUSIONS

This paper presents the most important results of a dynamic growth modeling study involving European Beech (Fagus sylvatica) forests, based on an extensive dataset provided by the Swiss Federal Institute for Forest, Snow and Landscape Research WSL, in Birmensdorf/Switzerland. The dataset includes 143 research plots, which had been re-measured between one and seventeen times, providing up to 16 interval measurements per research plot. In this model, the initial stand conditions at any point in time are defined by the three state variables dominant height, basal area and number of trees per hectare. According to this basic structure, the whole-stand model requires five stand-level inputs for simulation: the stand ages at the beginning and at the end of the projection interval as well as the initial dominant height, number of trees per hectare and stand basal area.

The model uses transition functions to project the state variables at any particular time. Two of these transition functions were derived using the generalized algebraic difference approach (GADA) and the third one was derived with the algebraic difference approach (ADA). These techniques ensure that base-age and path invariance properties provide consistent predictions. All the functions were fitted simultaneously using iterative seemingly unrelated regression and a base-age-invariant method that accounts for site-specific and global effects (Cieszewski et al., 2000).

The global whole-stand growth model proved to be robust for medium term projections of basal area. Considering the required accuracy in forest growth modelling, we can state that, based on the critical errors and relative RMSE, the model provides satisfactory predictions of basal area for time intervals up to 20 years.

Several single tree growth simulators have been tested during the past years for beech stands in Switzerland. However, their applicability depends on detailed stand structure information, which is not always available (Schmid et al., 2006). The whole-stand model presented in this study only requires five stand-level inputs for simulation: the age of the stand at the beginning and at the end of the projection interval as well as the initial dominant height, number of trees per hectare and stand basal area.

This model bridges the gap between the traditional yield tables and the desitable single tree models. It is preferable to yield tables in evaluating alternative silvicultural strategies in even-aged beech forests with varying densities and species combinations. One particular advantage is the fact that the weight and the type of a thinning may be specified quite easily using, for example, the basal area and stem number proportions removed in a thinning (refer to Gadow, 2006, p. 133). Moreover, the relatively simple structure of the growth model makes it suitable for decision support systems that enable forest managers to generate optimal management strategies and developing a reasonable design of a forested landscape in space and time (Gadow and Pukkala, 2008).

Acknowledgments

We are grateful to be able to use the extensive data base provided by the Swiss Federal Institute for Forest, Snow and Landscape Research.

References

All Tables

All Figures

	Figure 1 Map of all beech growth and yield plots (black dots = beech plots with tree height measurements, white dots = plots without height measurements). Map reproduced with permission of swisstopo (JA082265).
In the text

	Figure 2 H-Residuals versus H-Lag1-Residuals for model E1 (first row), G-Residuals versus G-Lag1-Residuals for model E1 (second row) and N-Residuals versus N-Lag1-Residuals for model E5 (third row) fitted using a first-order continuous time autoregressive error structure.
In the text

Figure 3 Upper left: dominant height curves for site indices of 18, 22, 26, 30 and 34 m at a reference age of 80 years for the GADA equation E1. Upper right: reduction in tree number curves for 400, 300, 2 200, 3 100, 4 000 trees per ha at a reference age of 40 years for the model (E5). Down: basal area growth rates for thinned (left) and unthinned stands (right) as affected by initial stem number and age overlaid on the trajectories of basal area.

In the text

	Figure 4 Plots of RMSE in percent (up) and critical error in percent (down) for basal area predictions against different prediction intervals. The rectangle indicates the maximum time interval for a required accuracy of 15% in basal area estimations.
In the text