Modeling dominant height growth of maritime pine in Portugal using GADA methodology with parameters depending on soil and climate variables

Abstract

• Background

Maritime pine (Pinus pinaster Ait.) is a common conifer species in Portugal that contributes significantly to the national economy. Accurate classification of forest productivity based on site index and height growth dynamics is the main basis for sustainable forest management of this species.

• Objectives

The main objective of this study was to develop a new dynamic site-dependent height–age model for the maritime pine in Portugal, using the generalized algebraic difference approach (GADA) methodology, and to explore possible improvements of the model´s performance by expanding its parameters as sub-functions of soil and climate variables.

• Methods

We tested for this purpose several dynamic equations, including anamorphic, polymorphic with common asymptote, and polymorphic with multiple asymptotes equations. The candidate models were fitted to a large set of stem analysis data, and tested on independent data from permanent sample plots.

• Results

The two best models with multiple asymptotes, one anamorphic and one polymorphic, showed similar performance; however, upon expanding the parameters as sub-functions of the climate and soil variables, the polymorphic model outperformed the anamorphic model, as well as other models previously used for the management of this species in Portugal. The results of this study also demonstrated that the maritime pine model, developed with stem analysis data, can accurately predict the dominant height growth measured on permanent sample plot data.

1 Introduction

The most common conifer species in Portugal is the maritime pine (Pinus pinaster Ait.). The maritime pine occupies 926 × 10³ ha (Tomé et al. 2007), which represents 27.2% of the total forested area in the country. Maritime pine has been present in Portugal for over 33,000 years (Figueiral 1995) and is present along the Atlantic coast, from the basin of the Tagus River, near Lisbon, to the northern Portuguese border. The species spreads inland to northern and central regions, ascending to altitudes of 700 to 900 meters where the effects of the Atlantic Ocean are prevalent (Oliveira et al. 2000). The area of maritime pine in Portugal increased from just over 250,000 ha to 1,300,000 ha in the period from 1902 to 1963 (Mendes 2007). Maritime pine has been an excellent choice for forestry restoration projects, including dune fixation in coastal areas and the restoration of forest cover in deforested areas, due to its rusticity, fast growth rate and ability to grow in uncolonized areas. As a result of forest fires, the area of maritime pine decreased to less than 1 million ha in the last 2 decades. In Portugal, maritime pine contributes significantly to the national economy, and is relevant in many industries. The species is commonly managed in rotations of 40–45 years (70–80 years in Leiria National Forest). Maritime pine is used for wood and resin production, as well as environmental and social services such as recreation and environmental protection.

Sustainable forest management requires a reliable estimation of wood quantity and quality. Growth and yield models are useful tools for these estimations and usually include a dominant height growth sub-model. For the range of stand densities usually used in forest management, dominant height growth is independent of stand density and is not affected by thinning; thus, dominant height is an important variable in characterizing the development of stands and estimating the growth potential of a site (Gadow and Hui 1999). The growth potential of a given site (site quality) is measured by the dominant height at a given age, commonly referred to as the site index.

Clutter et al. (1983) described three methods of constructing site index curves: the guide curve method, the difference equation method and the parameter prediction method. The guide curve method, a population approach based on average values (Zeger et al. 1988), has been used to develop site index curves when growth data were not available. The difference equation method and parameter prediction method are subject-specific approaches that recognize subject-specific variability and require re-measurement data from permanent plots or stem analysis data. Currently, dominant height growth over time is typically modeled with self-referencing functions (Northway 1985). These models predict dominant height at a time, t, as a function of height at an initial time, t _0. Parameters of self-referencing functions are estimated by simultaneously fitting all growth curves, without omitting individual site classes (Cieszewski et al. 2007). Among these three-dimensional site–height–age models, dynamic models, or base-age invariant equations, are more likely to provide numerous desirable properties that adequately describe height development (e.g., Cieszewski and Bailey 2000; Cieszewski 2003) such as curves through the origin, sigmoidal curves, polymorphism, variable asymptotes, equality of site index and predicted height at base-age and theoretical interpretability. The idea of base-age invariance property was formalized by Bailey and Clutter (1974), who introduced the algebraic difference approach (ADA). The ADA is implemented by assigning one site-specific (local) parameter in a base equation while all the others are global (common) parameters. Models based on this technique are either anamorphic or have a common asymptote, which is the main limitation of ADA (Cieszewski 2001). Cieszewski and Bailey (2000) generalized ADA and introduced the generalized algebraic difference approach (GADA), which allows for more than one site-specific parameter, usually two in practice. This approach allows the derivation of more flexible dynamic models, and is considered the most advanced method for modeling both polymorphism and multiple asymptotes.

Oliveira (1985) and Páscoa (1987) used the Schumacher (1939) equation along with the guide curve method to model dominant height growth of maritime pine in Portugal. Oliveira (1985) used a small data set of temporary plots to develop his site index curves, while curves developed by Páscoa were restricted to the Leiria National Forest. Tomé (2001) used a larger data set of permanent plots to fit an ADA equation based in the Lundqvist–Korf (Lundqvist 1957) function, where the parameter related to growth-rate was considered site-specific. Therefore, Tomé (2001) site index curves are polymorphic with a common asymptote.

The main purpose of this research was to develop a new well-behaved flexible base-age invariant dynamic equation fitted to a large set of stem analysis data for the maritime pine in Portugal. Furthermore, to gain insight on the development of site index curves, we tested the following hypothesis:

• 1. The predictive ability of a dynamic model of dominant height growth can be improved if some of the model parameters are expressed as a function of climate and soil variables.

• 2. Stem analysis data (individual tree data) can be used to develop site index curves for the prediction of dominant height growth (stand data).

The model should achieve the desirable properties describing height–age development in the sense of Cieszewski and Bailey (2000). Furthermore, the model should be applicable to the whole area of distribution of maritime pine in Portugal.

2 Data

Typically, data from the re-measurement of sample plots or from stem analysis are desirable for the development of site index curves (Cieszewski 2004). We used a stem analysis database from 247 dominant trees for model fitting. This database covered the entire range of maritime pine in Portugal. Trees were sampled from 500 m² temporary or permanent plots by selecting two or five dominant trees from each plot. If two trees were sampled, the method proposed by Tennent and Burkhart (1981) was used for tree selection. Tree ages ranged from 11 to 69 years old, and heights ranged from 6.7 to 29.8 m. Trees were felled, leaving stumps with an average height of 0.10 m. Felled trees were sectioned at the dbh level and at intervals of 2.0 m thereafter. Trees are usually not sectioned at the whorls level; thus, the true total height of a tree at an age based on the ring count at a crosscut is almost always underestimated (Carmean 1972). To correct for this bias, we employed the method of Carmean (1972) and the modification proposed by Newberry (1991). A total of 9,094 height/age pairs were used to fit the dynamic equations.

For model validation, we used an independent data set with 145 growth series obtained from re-measured permanent sample plots and thinning trials. Trials were located in sites from across the country, covering a broad range of maritime pine growth conditions. These trials were established between 1981 and 1998 in pure even-aged stands dating back to the 1960s and the 1970s. The plot size ranged between 500 m² and 1,000 m², and the number of re-measurements varied between two (1994 and 1999) and ten (yearly, from 1986 to 1992, and then in 1995, 1998 and 1999). The stem analysis and permanent sample plot data sets are summarized in Table 1.

Table 1 Data summary for fitting and validation data sets

Climate (averages from 1961 to 1990) and soil data were obtained from a GIS database with a resolution of 2 × 2 km². The data was based on digital maps from the Portuguese Agency for the Environment, and covered the entire country ( www.iambiente.pt/atlas/est/index.jsp ). Additionally, more detailed information (Daveau 1985) on the climate and soil of the sampling sites were obtained. We processed the information with ArcGis 9.2 from ESRI (ESRI Inc. 2006), where the data in the grid was spatially linked to the shape containing the coordinates of sampling locations. Table 2 presents descriptive statistics related to climate and soil characteristics of the sampling sites.

Table 2 Climate and soil variables used to expand the parameters of dynamic equations (values of climatic variables are annual mean values)

3 Methods

3.1 Approach for dominant height modeling

In this study, we used the generalized algebraic difference approach (GADA) to model the dominant height (hdom) of maritime pine. GADA allows several parameters to vary between sites (site-specific) as long as, by some algebraic transformation, they all can be expressed as a function of fixed or global parameters and just one varying parameter (Cieszewski 2003). The method consists of (Cieszewski et al. 2007): i) choosing a base equation to model the variable of interest (hdom in this case), ii) deciding which parameters of the function to relate with a theoretical site quality measure (X) and expressing this relationship through a mathematical equation, iii) solving the equation for X, and iv) obtaining the dynamic three-dimensional site–height–age model (Cieszewski et al. 2007) in the form of hdom = f (t, t ₀, hdom ₀) by substituting the solution of X in the explicit three-dimensional equation hdom = f (t, X) for the initial conditions, t₀ and hdom ₀. 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, GADA is equivalent to ADA, the first method used to develop self-referencing functions (Bailey and Clutter 1974). ADA has been successfully used by several authors (e.g., Tomé 1989; Elfving and Kiviste 1997; Amaro et al. 1998). Álvarez-González et al. (2010) recently used ADA for modeling tree mortality and GADA for modeling basal area and dominant height growth.

A total of 24 base-age invariant dynamic equations derived with the GADA approach (Cieszewski and Bailey 2000), were selected from the literature and fitted to the stem analysis data set (Table 3): four anamorphic equations (T1, T4, T8 and C10), six common asymptote polymorphic equations (T2, T3, T9, T6, T7 and M1) and 14 multiple asymptotes polymorphic equations.

Table 3 Base models and tested dynamic equations derived from the GADA approach (h = hdom; h ₀   = hdom ₀ )

3.2 Parameter estimation

An effective technique for parameter estimation in self-referencing functions should account for individual trends in the data (Cieszewski et al. 2000). In this study, we used the varying parameter (VP) method (also known as the dummy variable approach). With this method, all the data points are used to produce residuals, and arbitrary choices regarding measurement intervals are unnecessary. The procedure is considered a base-age invariant method, and produces unbiased estimates for base-age invariant equations. We implemented the VP method as in Cieszewski et al. (2000). In this technique, global and site-specific parameters are simultaneously estimated, resulting in the same parameter estimates regardless of the selected base age. The equations were fitted under the following form:

$$ hdo{m_{{ij}}} = f\left( {{t_{{ij}}},\,hdo{m_0},\,{t_0},\,\beta } \right) + {\varepsilon_{{ij}}} $$

(1)

In the fitting procedure, the variable hdom ₀ is substituted by $\sum\limits_{{i = 1}}^n {{S_i}\,{d_i}}$, a sum of terms containing a site-specific parameter, S _i , and a dummy variable for each tree, d _i , which is equal to 1 for the i ^th tree and 0 otherwise. Any initial height, corresponding to an arbitrarily selected initial age, t ₀ (t ₀  = 0 is not allowed), is used as starting value for S _i . The sum of terms collapses into a single parameter unique to each tree. With respect to the other terms in equation (1), hdom _ij is the dominant height of i ^th tree at an age t _ij , β is a vector of global parameters, and ε _ij is the model error.

Models must be developed with proper statistical procedures, and should not violate underlying assumptions (e.g., independent and identically distributed [iid] errors [ε _ij ] with a mean of zero, homogeneous variance and a normal distribution). Due to the longitudinal nature of stem analysis data used in model fitting, correlated errors within the same individual were expected. Therefore, we modeled the error autocorrelation, when present, using a stationary autoregressive structure of order p (AR(p) suitable for equally spaced measurements in time. With the AR(p) structure, the error term is expanded as follows:

$$ {\varepsilon_{{ij}}} = \varphi {_1}\,{\varepsilon_{{ij - 1}}}\,\,{d_1} + ... + \varphi {_p}{\varepsilon_{{ij - p}}}\,{d_p} + {u_{{ij}}} $$

(2)

where ε _ij is the error related to the j ^th measurement on the i ^th tree, φ _k are the k parameters to be estimated for the autoregressive process of order p, d _k is a dummy variable equal to 1 when j > k and equal to zero when j < k, and u _ij is white random noise.

We used the graphical representations of studentized residuals (rstd) versus lagged rstd to evaluate the order of the autoregressive structure. These residuals were obtained by initially fitting the models using PROC NLIN in SAS/STAT (SAS Institute, Inc. 2004a) without expanding the error term.

The lack of normality for errors observed in residual analysis (QQ-probability plots) was resolved by applying the Huber function (Myers 1986) for robust estimation. The VP method, including the corrections for autocorrelation of the residuals and for non-normality, was implemented in SAS/ETS using PROC MODEL (SAS Institute Inc. 2004b). We examined possible violations of the assumption of homoscedasticity by plotting standardized residuals versus predicted values of dominant height.

3.3 Evaluation of the dynamic equations

To compare the performance of the candidate dynamic equations and to select the best function for modeling dominant height growth of maritime pine, statistical and graphical analyses were conducted. Model evaluation was conducted in two phases:

The fitting performance of the models was evaluated by a statistic equivalent to the adjusted R-square in linear regression (MEF _adj ), the root mean square error (RMSE), and the bias (ē).

$$ ME{F_{{adj}}} = 1 - \frac{{n - 1\,\sum\limits_{{i = 1}}^n {{{\left( {{y_i} - {{\hat{y}}_i}} \right)}^2}} }}{{n - p\,{{\sum\limits_{{i = 1}}^n {\left( {{y_i} - \bar{y}} \right)} }^2}}} $$

(3)

$$ RMSE = \sqrt {{\frac{{\sum\limits_{{i = 1}}^n {{{\left( {{y_i} - {{\hat{y}}_i}} \right)}^2}} }}{{n - p}}}} $$

(4)

$$ \bar{e} = \frac{1}{n}\sum\limits_{{i = 1}}^n {\left( {{y_i} - {{\hat{y}}_i}} \right)} $$

(5)

where p is the number of parameters in the equation, n is the number of observations, y _i is an observed dominant height, and ŷ _i is the corresponding predicted dominant height.

The quality of fit does not necessarily reflect the quality of prediction, and the models should be tested with independent data (e.g., Myers 1986; Kozak and Kozak 2003). The dominant height growth model was developed with stem analysis data (tree data) but will be operationally applied to forest inventory data of dominant height (stand level data). It is therefore important to predict the behaviour of the model when operationally used. In this study, we used independent data from permanent sample plots for model validation. The data set was organized in the form of all possible intervals. To compare predicted and observed height, the model efficiency (MEF) was used:

$$ MEF = 1 - \frac{{\sum {{{\left( {{y_i} - {{\hat{y}}_i}} \right)}^2}} }}{{\sum {{{\left( {{y_i} - \bar{y}} \right)}^2}} }} $$

(6)

where y _i and ŷ _i are observed and estimated dominant height values respectively.

To characterize model error, we used two simple recommended criteria (Soares et al. 1995) that provide a summary of the overall model performance: the average model bias $\left( {\sum {\left( {{y_i} - {{\hat{y}}_i}} \right)} \,/\,n} \right)$ and the mean absolute difference $\left( {\sum {\,\left| {{y_i} - {{\hat{y}}_i}} \right|} \,/\,n} \right)$, where n is the total number of observations in the validation data set. The average model bias measures the error when several observations are combined by totaling or averaging, and the mean absolute difference measures the average error associated with a single prediction. Graphical representations of these statistics as versus projection length interval were also used as validation criteria. We analyzed the biological consistency of the models by observing the magnitude and sign of the parameter estimates, as well as the estimated upper limit value of dominant height for an age of 100 years ($h\widehat{d}o{m_{{100}}}$). At this age, maritime pine has reached its maximum height. Thus, satisfactory model performance is indicated by an accurate prediction of $h\widehat{d}o{m_{{100}}}$, even if the model has a higher upper asymptote. Although maritime pine can live up to 200 years old, only a few isolated specimens actually reach this age. In practice, most trees live for approximately 80–100 years (Lanier 1986). We used the highest value of hdom (hdom ₀   = 29.8; t ₀   = 59) observed in the fitting data set, corresponding to the highest site index value (SI ₅₀   = 28 m) estimated according to the equation of Tomé (2001), to compute $h\widehat{d}o{m_{{100}}}$.

We compared the predictive ability of the best models to existing Portuguese site index models, including the equation of Tomé (2001). This equation was developed with the data set that was used in this study for validation purposes. Thus, this model can be used as a standard for comparison purposes with other models because its predictive performance is high in the data set to which it was fitted.

3.4 Expanding model parameters with climate and soil data

Soil and climate gradients are among the most important factors that explain the growth capacity of forest species (Oliver and Larson 1996). We studied possible improvements in dominant height estimation by relating parameter estimates to climate and soil variables, as Wang et al. (2007) and Bravo-Oviedo et al. (2008) did. The best models were fitted with parameters expressed as combinations of several variables (Table 2). In each combination, one parameter was related to soil and the other was associated with climate variables. An ordinal scale was used for variables summer type (SUMMER) and winter type (WINTER). Dummy variables were employed for soil type (ST) and soil parent material (SM) (see Table 2).

4 Results

Analysis of the plots of studentized residuals (rstd) versus lagged rstd showed that the first-order autoregressive structure (AR(1)) was appropriate for modeling the error term.

Table 4 shows goodness-of-fit statistics resulting from model fitting. For simplification, only the results of the ten best models are presented.

Table 4 Fitting and validation statistics for the ten best models

Among the models that best fit the stem analysis data set, one equation was anamorphic (T4) and three were common asymptote polymorphic equations (T5, T6 and T7) from Tomé (1989). Furthermore, six multiple asymptotes polymorphic equations (C4 to C7, C14 and C15) from Cieszewski (2004) and Cieszewski et al. (2006) were also included. Generally, differences in global fitting results were minor.

The statistics computed in the validation phase, using independent data from permanent sample plots, are also presented in Table 4. Again, only the ten best models are shown. Seven of these models are common to the group of models on the left side of Table 4, namely the anamorphic equation, T4, the common asymptote polymorphic equation, T5, and the multiple asymptotes polymorphic equations C6, C7, C14 and C15. Three models provided accurate predictions of dominant height, although they did not provide the best fit. Specifically, these models include two anamorphic equations (T1 and C10) from Tomé (1989) and Cieszewski (2002) respectively, and one multiple asymptotes polymorphic equation (C3, Cieszewski, 2004).

Only a few rare specimens of maritime pine in Portugal have been reported to achieve a total height greater than 35 m. The Portuguese Forest Service reports a 200-year-old isolated specimen of Pinus pinaster in Leiria National Forest with a height of 40 m and a perimeter of 4.51 m at the dbh level. The maximum height estimated at an age of 100 years was used as a reference to evaluate asymptotic behavior of the models. For this estimate, a value between 35 to 40 m was considered adequate, even if the model presented a maximum asymptote above this point. Six equations in Table 4 predicted an $h\widehat{d}o{m_{{100}}}$ > 35 m. However, equation T7 was excluded from this set because the predicted value of $h\widehat{d}o{m_{{100}}}$ was near the lower limit of the acceptable range. Thus, five equations were selected for further analysis. Parameter estimates for these equations are presented in Table 5. All parameters were significantly different from zero (p < 0.0001), and displayed consistent magnitudes and appropriate signs to model biological growth.

Table 5 Parameter estimates for the five best equations modeling hdom growth

Figure 1 shows the bias and precision of models projecting dominant height. Models T4 and C5 were the most accurate, presenting similar bias and precision at different intervals. Values of the average bias and the mean absolute difference confirm the better prediction ability of these models (Table 4). These two equations were derived from the Lundqvist–Korf (Lundqvist 1957) function by applying the GADA approach. Despite the fact that T4 is anamorphic and C5 is polymorphic with multiple asymptotes, both presented an equal value of model efficiency (98.2%) and $h\widehat{d}o{m_{{100}}}$. Site index curves obtained with these two models were similar (Fig. 2).

Fig. 1

Bias and precision of dominant height time projections

Fig. 2

Site index curves from C5 and T4 for a base-age of 50 years

In both models, the potential improvement in dominant height prediction was analyzed by expanding the parameters as a function of climate and soil characteristics.

The parameters b ₁ and c of the dynamic equation C5 (Table 3) were related to soil and climate variables. With respect to soil type, i-1 indicator variables (ST_i, i = 1, ..., 5) were used to expand the parameter c (i is the i ^th soil type). The two best combinations were:

$$ {b_1} = {b_0} + {b_2}\,\frac{{INS}}{{ELEV + 1}}; c = {c_0} + \sum\limits_{{i = 1}}^4 {{c_i}\,S{T_i}} $$

(7)

$$ {b_1} = {b_0} + {b_2}\,\frac{{P\,\sqrt {T} }}{{WINTER}}; c = {c_0} + \sum\limits_{{i = 1}}^4 {{c_i}\,S{T_i}} $$

(8)

Upon fitting these models, it was observed that rankers (ST ₂ ) and calcic cambisols (ST ₃ ) could be combined into one group. Orthic podzols (ST ₄ ) and dystric cambisols (ST ₅ ) were combined with humic cambisols (ST ₁ ) to form another group, which was used as the base level. In models C5_B and T4_B, humic cambisols (ST ₁ ) were detached from the base group and considered as a third group.

Expanding b and c of the dynamic equation T4 (Table 3) with soil and climate variables, models T4_A and T4_B were obtained. Model efficiency, bias and precision for C5, T4 and the corresponding expanded models are presented in Table 6. Table 7 shows the values of parameters of the expanded versions of T4 and C5. The predictive ability of the models was improved due to the introduction of these variables into the models parameters.

Table 6 The predictive ability of C5 and T4, expanded models and other relevant dominant height growth equations used in Portugal

Table 7 Parameter estimates for the expanded versions of C5 and T4

5 Discussion

This study evaluated several dynamic equations for modeling dominant height growth of maritime pine in Portugal. Dynamic site equations, in contrast to static site equations, directly estimate height and site index from any other height and age, being base-age invariant (Cieszewski 2001). Because dynamic equations can predict both site index and dominant height, they are more parsimonious and flexible than static equations, and are able to define broader ranges of curve shapes. As stated by Cieszewski (2004), the search for methods to increase flexibility in site index equations has been an important goal of forest biometricians (e.g., Borders et al. 1984). The generalization of the algebraic difference approach introduced by Cieszewski and Bailey (2000) was considered a major step forward in dynamic equation modeling (Cieszewski et al. 2007) because these methods increased flexibility, and allowed for the development of equations that account for both polymorphism and multiple asymptotes, which are considered the most advanced form of dynamic equations. In this study, anamorphic and polymorphic equations with common asymptote and multiple asymptotes were compared. However, not all equations that performed best in the fitting and validation phases were polymorphic equations with multiple asymptotes. For example, the anamorphic equation T4 fitted the data well and displayed similar predictive ability to the best multiple asymptote polymorphic model C5. As shown in Table 4, other anamorphic and common asymptote polymorphic equations performed well. Many authors have reported the superiority of multiple asymptotes polymorphic equations to describe the growth patterns in their data in comparison to either anamorphic or common asymptote polymorphic equations (e.g., Cieszewski 2003; Dieguez-Aranda et al. 2005; Cieszewski et al. 2006, 2007). However, this superiority was not so evident in our study. The two best models, C5 and T4, displayed similar accuracy in dominant height estimation. Surprisingly, the performances of T4 and C5 were similar, indicating that polymorphism in the height growth of maritime pine in Portugal is not clear. Developing methods that can efficiently identify anamorphic or polymorphic trends in site index data is an important area of future research. Nevertheless, in polymorphic model C5, two parameters are considered to vary with the site conditions, which is an advantage over the T4 model. Thus, C5 accounts for concurrent polymorphism and multiple asymptotes, two important properties frequently required from site index curves. Furthermore, in GADA methodology, X may incorporate variables related to climate, soil, photosynthesis and related information, as well as genetic parameters (Cieszewski and Bailey 2000). Thus, another advantage of C5 over T4 is the greater flexibility in introducing information regarding the tree’s environment, allowing more site-specific parameters. Furthermore, equation C5 is also biologically more realistic than T4, which is often used as a criterion for selecting models that present similar statistical performance. When parameters of equations C5 and T4 were expanded as sub-functions of soil and climate variables, the accuracy of height estimations further increased. In version A of the models, parameter c was considered to be a linear function of soil type (ST), and the growth-rate related parameter (b ₁ in C5_A and b in T4_A) was considered a linear function of the ratio between insolation (INS) and elevation (ELEV). Insolation is a measure of solar radiation that reaches a given surface area over a specific time period. In this study, average annual values of insolation from 1961 to 1990 were used. Solar radiation energy is an important factor for tree growth, being the primary source of energy for the synthesis of organic material (Larcher, 2007). However, at higher elevations, the climate becomes progressively cooler and moister (Oliver and Larson 1996). Therefore, models such as C5_A and T4_A are physiologically and ecologically sound, because their parameters depend on climate and soil characteristics. In the B version of the models, parameter c was related to soil type, and the growth-rate parameter was directly related to precipitation (P) and temperature (T), and inversely related to the winter pattern (WINTER). B versions of the models were less precise (only slightly) than the A versions, but displayed identical model efficiency and lower bias. Oliveira et al. (2000) considered a mean annual precipitation of 800 mm (with a minimum of 100 mm during the annual summer dry period) and mean annual temperature values of 13 to 15°C to be ideal conditions for maritime pine growth. Furthermore, the species is sensitive to long periods of frost. Thus, these climate variables were included in models C5_B and T4_B, where characteristics of winter are related to frost patterns. Bravo-Oviedo et al. (2008) reported that the relation of parameters from a previously selected dynamic equation (Bravo-Oviedo et al. 2007) to climatic and soil characteristics led to an improvement in site index modeling of Mediterranean maritime pine in Spain. The authors found that a model incorporating additive effects of seasonal precipitation and temperature ($PR\sqrt {T}$), as well as dry summer period length and soil type (dolomite versus other soil types) performed particularly well. Application of this climate-based model led to more accurate site index predictions at a regional level.

A comparison of the predictive ability of the expanded versions of C5 and T4 revealed minimal differences, which was similar to the results of the unexpanded models. However, the validation statistics in Table 6 show that the expanded versions of C5 have a slight advantage (C5_A over T4_A and C5_B over T4_B), most probably due to greater flexibility. Moreover, expanding T4 with soil and climate variables caused the model to lose its anamorphic properties.

Bias and precision of the two versions of C5 and T4 are presented in Fig. 3 by age and site index classes. Model C5_B was slightly less biased than C5_A in almost all site index classes and age classes (15 to 25 years). The precision is slightly higher in C5_A than in C5_B, but differences were minimal. Models T4_A and C5_A performed similarly. However, C5_A displayed a slight advantage. Model T4_B presented the poorest performance. All four models outperformed other equations used in the modeling of dominant height of maritime pine in Portugal (Table 6), including the Tomé (2001) equation, which was taken as a benchmark for the performance of other models. Models that include parameters expanded with sub-functions of climate and soil variables are more robust than non-expanded ones to geographical variability induced by climate and or edaphic conditions, which can cause type II polymorphism in the sense of Krumland and Eng (2005).

Fig. 3

Bias and precision of C5_A and C5_B by age and site index class. * Site index (SI) computed with the Tomé (2001) equation

Some authors claim that if empirical models are provided with climate variables in a consistent way, they can play an important role in assessing the effects of climate change on the growth trajectories, helping to take adaptive measures in forest management. This assumption should be used with caution, as it may lead to erroneous conclusions. In this study, the climate and soil variables included in the expanded models aim to localize the growth curves at the regional level, and not to predict the response to climate change. Empirical models are specific to the sites where the measurements made to determine parameter values have been taken, and should not be applied to different regions and conditions (Landsberg, 2003). Particularly, the fitted models include just a few climatic or climate-related variables expressed as long-term average climatic data that, even if combined with soil type, do not make the models sensitive to inter-annual or seasonal variability, nor are able to mimic the complex interconnectedness of the effects of site variables. Moreover, no site condition thresholds that could result in irreversible impacts in tree growth are considered. Prediction of tree dominant height responses to climatic change requires an understanding of the relationships linking site conditions to tree growth that is not present in the proposed models. Taking these limitations into account, the proposed models might be used to explore, only in the short term, the impact of inter-annual climatic variations.

Due to the nature of the climatic variables included in the models, equation C5_B may better accommodate short-term changes in climate conditions than C5_A or T4_A. Furthermore, C5_B is more useful for practical applications. Based on the above analysis, model C5_B is proposed as the best model for the construction of a new set of site index curves for maritime pine in Portugal.

6 Conclusions

The objective of this study was to develop a set of site index curves for the maritime pine in Portugal based on a large stem analysis data set. Testing if expanding the model parameters as a function of soil and climate variables would improve model performance, and also if models developed with stem analysis perform well when applied to permanent plot data, were also objectives of the study.

The selected model, a multiple asymptotes polymorphic equation from Cieszewski (2004), modified to include parameters relating to climate and soil characteristics, is given by the following expressions:

$$ hdom = \exp \,\left( {{X_0}} \right)\,\exp \,\left( { - \left( {{b_1} + \left( {1\,/\,{X_0}} \right)\,} \right)\,{t^{{ - c}}}} \right) $$

(9)

$$ {X_0} = 0.5\,\left( {{b_1}\,\,t_0^{{ - c}} + \ln {H_0} + {F_0}} \right) $$

(10)

$$ {F_0} = {\left( {{{\left( {{b_1}\,t_0^{{ - c}} + \ln {H_0}} \right)}^2} + 4\,\,t_0^{{ - c}}} \right)^{{1/2}}} $$

(11)

where ${b_1} = 9.2562 - 0.0006\,\frac{{P\,\sqrt {T} }}{{WINTER}}$

and $c = 0.5394 - 0.0918\,S{T_1} - 0.1767\,\left( {S{T_2} + S{T_3}} \right)$ (ST ₁ , ST ₂ , and ST ₃ are dummy variables for humic cambisols, rankers and calcic cambisols, respectively).

This study indicates that the performance of dominant height growth models can be improved by expanding parameters to include soil and climate characteristics. The proposed model, polymorphic and with parameters depending on soil and climate variables, improved site index modeling of maritime pine in Portugal and is consistent with ecological and physiological concepts governing the tree growth. Furthermore, the model can accommodate short-term changes in climatic conditions. The predictive ability of dominant height growth models developed with stem analysis data is often questioned when operationally applied to plot data. This study reveals that models developed with stem analysis data can perform well, and are comparable to models developed with permanent plot data.