Developing a dynamic growth model for maritime pine in Asturias (NW Spain): comparison with nearby regions

Abstract

∙Key message

A dynamic growth model was developed for maritime pine in Asturias. During the evaluation process, a stand volume ratio function proved the best of two alternative methods for estimating merchantable volume. Comparison of the developed model with existing models for nearby regions showed that a single model may suffice for the whole of the NW Iberian Peninsula.

∙Context

Maritime pine is one of the most important tree species in NW Spain. There was no existing dynamic growth model for this species in Asturias.

∙Aims

To develop a dynamic growth model for maritime pine in Asturias, by evaluating two different methods of estimating volume (a disaggregation system and a stand volume ratio function), and to compare the developed model with existing models for Galicia and northern Portugal are the goals of this study.

∙Methods

The dynamic model is based on the state-space approach, in which three state variables characterize the stand at any point in time: dominant height, number of stems per hectare and stand basal area. The transition function for the first variable was developed on the basis of stem analysis data in a previous study, while the corresponding functions for the last two variables were simultaneously fitted with data obtained from successive measurements of permanent plots. An appendix outlining the implementation of a stand growth simulator in the R environment is included to facilitate model use and evaluation.

∙Results

When the whole model was used to project the stand conditions, the stand volume ratio function performed best, yielding a root mean square error of 22.4 m³ ha⁻¹ and a critical error of 18.4 %. Comparison with models developed for other regions revealed both similarities and differences, some of which may be attributed to an unequal distribution of the available data in age and site quality classes.

∙Conclusion

The proposed dynamic growth model provided accurate results, and comparison with other region-specific models showed that a single dynamic model may suffice for the whole of the NW Iberian Peninsula.

1 Introduction

The growing stock of maritime pine (Pinus pinaster Ait.) is the largest in Spain, currently representing 15 % of the timber volume and 27 % of the annual harvested volume (MAGRAMA 2010). The wood is commonly used to produce sawn timber and also in the pulp or wood-based panel industries, depending on log dimensions (Sanz et al. 2006). Maritime pine covers a total area of 22,500 ha in the region of Asturias (NW Spain) (5 %, MAGRAMA 2012).

Because of the importance of this species in Asturias, several studies have addressed the following aspects: site quality (Álvarez-Álvarez et al. 2011; Arias-Rodil et al. 2015a), edaphic factors and nutritional status (Afif Khouri et al. 2009), tree biomass (Canga Líbano et al. 2009) and diameter distribution (Gorgoso Varela et al. 2009). However, no growth and yield models have yet been developed for the species. Such models are very useful for forest practitioners as they enable stand growth simulation for different management practices. This type of information is essential for decision-making purposes in the context of sustainable forest management (Vanclay 1994, p. xv). Empirical models are widely used as practical tools in forest management, as they require few and easy-to-measure (inexpensive) input variables and provide accurate yield predictions. However, as they are not linked to the underlying growth processes, they do not take into account sudden disturbances (e.g. extreme weather events) or long-term environmental changes (i.e. climate change). A detailed review of model types can be found in Mäkelä et al. (2000), Monserud (2003) and Pretzsch et al. (2008).

Empirical growth models can be classified into three types (Davis et al. 2001, p. 185): (i) whole-stand models, (ii) size-class models and (iii) individual-tree models. The differences between these models are based on the level of resolution of predictions and data requirements. Model types (ii) and (iii) are useful as research tools (they describe stand behaviour), but over-parameterization is common and limits the accuracy and precision of predictions at forest management level (García 2003). In addition, collection of individual-tree and size-class data is more expensive than collection of data at stand level. In the special case of the management of single-species plantations, whole-stand models therefore represent a good compromise between accuracy and generality (García 2003). They use stand variables as inputs and outputs, although some of the models may include a mathematical disaggregation of the stand-level inputs to estimate the number of stems in different diameter classes (referred to as diameter distribution models by Burkhart and Tomé (2012), p. 234).

Volume prediction to any merchantable limit (for different timber assortments) is commonly achieved (Burkhart and Tomé 2012) by using (i) a stand volume ratio function, which provides the volume to a top diameter limit (Barrio-Anta et al. 2008) or to both a top diameter limit and a diameter class (Amateis et al. 1986), or (ii) a disaggregation system, based on the joint use of a diameter distribution, a height–diameter (h–d) relationship, and a stem taper function, to estimate the merchantable volume according to top diameter limits and log lengths. The latter approach has been used in many models (Diéguez-Aranda et al. 2006; Castedo-Dorado et al. 2007; Gómez-García et al. 2014a), and it is the approach most commonly used to estimate volume within a dynamic growth model. Nevertheless, there are several disadvantages associated with this method, such as possible computational inefficiency (it depends on several functions, some of which may involve iterative procedures) and the inclusion of several submodels, with their corresponding estimation errors. Therefore, it is useful to evaluate whether a stand volume equation can yield results that are as accurate as those produced by the disaggregation system, but in a more computationally efficient way.

Growth and yield models have already been developed for maritime pine plantations in nearby Atlantic regions: the models developed for the species in Galicia are summarized in Diéguez-Aranda et al. (2009) and for northern Portugal (Tâmega Valley, model termed ModisPinaster) in Fonseca (2004) and Fonseca et al. (2012). These models provide information about the behaviour of this species in each region, enabling comparison with behaviour of the species in Asturias. We consider such comparison of interest, as it represents another way of validating the proposed dynamic growth model.

The main objectives of this study were as follows: (i) to develop a dynamic growth model for maritime pine in Asturias and to evaluate its performance in total and merchantable stand volume estimation, (ii) to compare the prediction of total and merchantable stand volume by using both a stand volume ratio equation and a disaggregation system and (iii) to compare the performance of the proposed dynamic growth model with those of existing models for Galicia and northern Portugal, as regards projection of stand variables, estimation of total stand volume and determination of optimal biological rotation age (the age of maximum mean annual increment—MAI—of total stand volume).

2 Material and methods

2.1 Data

The data used to develop the model were obtained from two networks of plots installed in pure, even-aged maritime pine stands: (i) 74 permanent plots and (ii) 18 plots of a thinning trial. The permanent plots of the first data source were established and measured in 2007. The plots were installed throughout the area of distribution of the species in Asturias (mainly in the NW of the region), covering the existing range of ages, stand densities and site qualities. The plot size ranged from 700 to 900 m² to include a minimum of 30 trees. A second measurement was made on a subset of 58 of the 74 initially established plots (in 2011 and 2012), as some disappeared as a result of forest fires or clear-cutting.

The 18 plots corresponding to the second data source are located in 6 sites (3 plots of 1000 m² per site), in which each plot was subjected to a different thinning treatment (no thinning, control; light low thinning and heavy low thinning). Three measurements were carried out (in 2009, 2011 and 2013), which implies two available growth intervals for each plot.

In each plot corresponding to the first and second data sources, diameter at breast height (d, cm, at 1.3 m from the ground) and total height (h, m) were measured to the nearest 0.1 cm and 0.1 m, respectively, in all trees. Descriptive variables were also recorded for each tree, e.g. if they were alive or dead. The following stand variables were calculated: age (t, years), dominant height (H, m, defined as the mean height of the dominant trees, i.e. the 100 thickest trees per hectare), site index (S, defined as the dominant height of the stand at a reference age of 20 years), number of stems per hectare (N), stand basal area (G, m² ha⁻¹), and total (V, m³ ha⁻¹) and merchantable (V _i ) volumes to different top diameter limits (from 1 to 35 cm by 1-cm intervals). The stand volumes were obtained by aggregation of the corresponding tree volumes estimated using the stem taper function fitted by Arias-Rodil et al. (2015b). In addition, a dominant height projection function was developed by Álvarez-Álvarez et al. (2011) on the basis of stem analysis data from 146 trees located close to the plots of the first data source (see Álvarez-Álvarez et al. (2011) for the data description; observed and predicted trajectories in Fig. 1, top left), which was therefore used as the transition function for dominant height and for computing the site index of each plot-inventory combination. Its expression is shown in Eq. 1 and the corresponding parameters are shown in Table 2. Site index values were then averaged by plot, assuming that site quality is constant over time.

$$ H_{2} = \frac{b_{1}}{1 - \left( 1 - \frac{b_{1}}{H_{1}}\right) \left( \frac{t_{1}}{t_{2}}\right)^{b_{3}}} $$

(1)

where H ₂ is dominant height (m) at age t ₂ (years), obtained from dominant height H ₁ at age t ₁, and b ₁ and b ₃ are the parameters of the ADA formulation of Hossfeld (1822) model, presented by McDill and Amateis (1992).

Fig. 1

Observed trajectories (solid grey lines) of dominant height (H, left) and stand basal area (G, right) by region, and the corresponding height predictions (site indices of 7, 11, 15 and 19 m) and stand basal area (values of 10, 25, 40 and 55 m² ha⁻¹ at 20 years), provided by the region-specific model (solid black lines), the submodels for Asturias (dashed black lines), and the submodels for Portugal (dotted black lines). The Portuguese stand basal area submodel depends on stem density at both initial and projection age, as well as on stand basal area at the initial age; the graph was obtained from the mortality curve that passes through 1000 stems ha⁻¹ at 20 years (the mortality curve also depends on dominant height, which was obtained from the height growth curve passing through 11 m at 20 years)

Transition functions serve to describe the natural evolution of stands, and therefore growth intervals used for fitting should not include silvicultural treatments (e.g. thinning) or random disturbances (e.g. fire, wind damage). In accordance with this criterion, 37 growth intervals in the first data source and 8 in the second (4 between the first and the second inventories and 4 between the second and the third) were disregarded. Thus, the transition functions were finally fitted for 88 plot-inventory combinations (35, 39 and 14 plots for the first, second and third inventories, respectively). In addition, all the 186 plot-inventory combinations available (92, 76 and 18 for the first, second and third inventories, respectively) were used to develop part of the disaggregation system (diameter distribution and height–diameter relationship) and the stand volume ratio function. Summary statistics of the stand variables used in model development are given in Table 1.

Table 1 Data summary

2.2 Dynamic model

The dynamic whole-stand model developed in this study is based on the state-space approach (García 1994), which is suitable for systems that evolve over time, such as forest stands. The current state of a stand can be defined by a list of state variables, which can be projected to the future by the transition functions. The idea is to characterize the state of the system at any point in time so that the future state does not depend on the past state. Three state variables are often used (dominant height, number of stems per hectare and stand basal area, e.g. Diéguez-Aranda et al.2006; Castedo-Dorado et al.2007; Álvarez-González et al.2010; Gómez-García et al.2014a), although for high-intensity treatments (e.g. high-height pruning, heavy thinning), use of a fourth variable, such as a measure of stand closure, is recommended to account for site occupancy and to take into account the response to these treatments (e.g. García 2011, 2013, García et al. 2011). In this study, we considered single-species stands, moderately thinned from below, and under these conditions, H, N and G are able to provide a good description of the stand condition at any age. Three transition functions enable projection of the state variables over time. According to García (1994), transition functions must possess certain desirable properties: consistency (no change for zero elapsed time), path-invariance (the result of projecting from t ₁ to t ₂ and then from t ₂ to t ₃ must be the same as projecting directly from t ₁ to t ₃) and causality (a change in the state can only be influenced by inputs within the relevant time interval). Once the state variables are known for a given time, total and merchantable volumes can be estimated in different ways. One way is to use a stand volume ratio function (e.g. Tewari et al. 2014) to estimate total and merchantable volume directly from the state variables. An alternative approach is to use a disaggregation system (e.g. Diéguez-Aranda et al.2006, Castedo-Dorado et al. 2007, Gómez-García et al.2014a): diameter distribution models generate the number of stems in diameter classes, while generalized height–diameter models predict the height for the average tree of each class, and taper functions are used to compute and classify the volume according to top diameter limits and log lengths, which are specified by market requirements.

In the following sections, we explain how we developed the transition functions, the stand volume ratio function and the disaggregation system. As already mentioned, Álvarez-Álvarez et al. (2011) have already developed the height growth curves using stem analysis data (1), and therefore, these were not re-fitted in this study. Most of the submodels were fitted by the ordinary least squares (OLS) technique, by using the nls function of R (R Core Team 2015).

2.2.1 Mortality and stand basal area growth functions

The algebraic difference approach (ADA—Bailey and Clutter 1974) and its generalization (GADA—Cieszewski and Bailey2000) were used to develop the transition function for number of stems per hectare and stand basal area. The transition function for the decrease in number of stems considers only the natural mortality (i.e. that caused by competition for light, water and soil nutrients—Peet and Christensen1987). Preliminary analysis showed that the best model for describing mortality was that of Tomé et al. (1997):

$$ N_{2} = N_{1} \exp{(a_{0} (t_{2} - t_{1}))} $$

(2)

where N ₂ is the projected number of stems per hectare at age t ₂ (years), N ₁ and t ₁ are respectively the initial values of number of stems and age and a ₀ is the model parameter.

Use of the ADA and GADA approaches implies that for prediction, the initial value of the variable must be known at a given age. This value is commonly known for the number of stems (which are easily counted), but not for the stand basal area. If this variable is not known, an initialization equation is used to predict stand basal area from other stand variables. Therefore, we developed a stand basal area growth system comprising two functions: one for projection and another for initialization. For modelling stand basal area growth, previous analyses showed that the best model was the GADA formulation of Hossfeld (1822) model (3). For the initialization function, previous analyses showed that basal area is allometrically related to S, N and t (4). In addition, the model was fitted with data from plots younger than 20 years, given that its predictive capacity was found to be much lower when the whole data set was used. Thus, for stands older than 20 years, we assumed that stand basal area should be obtained from inventory data.

$$ G_{2} = \frac{X_{0}}{1 + \frac{b_{2}}{X_{0}} t_{2}^{-b_{3}}} $$

(3)

where $X_{0} = \frac {1}{2} \left (G_{1} + \sqrt {{G_{1}^{2}} + 4 b_{2} G_{1} t_{1}^{-b_{3}}}\right )$; G ₂ is the projected stand basal area (m² ha⁻¹) at age t ₂ (years); G ₁ and t ₁ the initial values of stand basal area and age, respectively; and b ₂ and b ₃ the model parameters (b ₁ was non-significant).

$$ G_{1} = a_{0} S^{a_{1}} N^{a_{2}} t^{a_{3}} $$

(4)

In the process of fitting both transition functions, we used the base-age invariant dummy variables method proposed by Cieszewski et al. (2000), which estimates site-specific effects under the assumption that data measurements always include measurement and environmental errors (both on the left- and right-hand sides of the model) that must be modelled. Previous fitting showed that the residuals of both the transition function for stand basal area and number of stems were correlated (Pearson’s moment product correlation of 0.24). To account for this correlation, we fitted the models simultaneously by using the seemingly unrelated regression (SUR—Zellner1962) technique for nonlinear models. This method is completed in two steps, which are based on the OLS technique: (i) separate fitting of each of the equations considered in the system (the two transition functions in this case) and (ii) re-fitting simultaneously all equations of the system considering the correlation between residuals obtained in the first step.

The first step of SUR fitting was accomplished by the OLS technique, while for the second step, the residual sum-of-squares computation of the system was implemented in a function and then minimized by the nlm function of R (R Core Team 2015).

2.2.2 Volume estimation

Stand volume ratio function

Predicting the merchantable volume (to a top diameter limit) as a ratio of total volume was originally proposed by Burkhart (1977) at tree level, for which a total volume and a ratio equation must be fitted. Barrio-Anta et al. (2008) applied this approach to the estimation of merchantable stand volume. As in the latter study, we considered a total stand volume equation based on an allometric relationship to stand basal area and dominant height, and a ratio equation depending on quadratic mean diameter and a top diameter limit. Both were combined in the same equation, given that total stand volume becomes a special case of the stand volume ratio equation when the top diameter limit is equal to zero (Gregoire and Schabenberger 1996). The final equation has the following form:

$$ V_{\text{i}} = V R_{\text{i}} = a_{0} G^{a_{1}} H^{a_{2}} \exp\left( b_{0} D_{\text{g}}^{b_{1}} d_{\text{i}}^{b_{2}}\right) $$

(5)

where V _i is the stand volume (m³ ha⁻¹) to a top diameter limit d _i (cm), G the stand basal area (m² ha⁻¹), H the dominant height (m), D _g the quadratic mean diameter (cm) and a ₀-a ₂ and b ₀-b ₂ are the model parameters. Note that $V = a_{0} G^{a_{1}} H^{a_{2}}$ represents the total stand volume equation and $R_{\text {i}} = \exp \left (b_{0} D_{\text {g}}^{b_{1}} d_{\text {i}}^{b_{2}}\right )$ the ratio equation.

Disaggregation system

After selecting an appropriate functional form to characterize the diameter distribution, we obtain the parameter estimates and predict the number of stems per diameter class. The height of the average tree in each diameter class is predicted using the h–d relationship. Total and merchantable tree volume as dependent on top diameter and log length are predicted using a stem taper function and aggregating the individual-tree results to the whole stand.

Diameter distribution

We used a two-parameter Weibull function to model the diameter distribution, given that it is the most frequently used in forest growth models because of its flexibility and simplicity (Maltamo 1995; Kangas and Maltamo 2000; Torres-Rojo et al. 2000). The expression of the two-parameter Weibull density function is as follows:

$$ f(x) = \left( \frac{c}{b}\right) \left( \frac{x}{b}\right)^{c-1} \exp\left( - \left( \frac{x}{b} \right)^{c} \right) $$

(6)

where x is the random variable, and b and c (both positive) are the scale and shape (skewness) parameters, respectively.

The parameters of the Weibull function can be obtained by different methods, which can be classified in two groups: parameter prediction and parameter recovery (Hyink1980; Vanclay1994 p. 23). Considering the scope of this study, parameter recovery was used as it proved the best method in several previous studies (Cao et al. 1982; Reynolds et al. 1988; Torres-Rojo et al. 2000). The parameters are directly obtained from percentiles (Cao and Burkhart 1984) or moments (Newby 1980; Burk and Newberry 1984) of the diameter distribution, estimated from its relation with stand variables. Recovering parameters from the moments (known as the moments method) is the only method that directly ensures 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 (Hyink 1980; Kangas and Maltamo 2000; Torres-Rojo et al. 2000); we therefore chose to use this method.

The function parameters can be recovered from the first two moments of the diameter distribution (mean and variance) (Cao et al. 1982). The arithmetic mean diameter (D _m, cm) corresponds to mean, and the variance (var) is estimated as var = ${D_{\text {g}}^{2}} - {D_{\text {m}}^{2}}$ (D _g, quadratic mean diameter, cm). The known values of D _g and D _m can be used to obtain var and Eq. 7 can be solved iteratively for parameter c. Finally, parameter b is computed directly from Eq. 8 using previously obtained values of D _m and the parameter c.

$$ \text{var} = \frac{{D_{\text{m}}^{2}}}{{\Gamma}^{2} (1 + 1/c)} \left[ {\Gamma} \left( 1 + \frac{2}{c}\right) - {\Gamma}^{2} \left( 1 + \frac{1}{c}\right) \right] $$

(7)

$$ b = \frac{D_{\text{m}}}{\Gamma (1 + 1/c)} $$

(8)

where D _m the arithmetic mean diameter (cm), b and c the Weibull parameters and Γ the Gamma function.

The quadratic mean diameter (D _g) can be obtained directly from N and G $\left (D_{\text {g}} = 100 \sqrt {\frac {4 G}{\pi N}}\right )$, but D _m remains unknown. However, it can be modelled through a relationship with the quadratic mean diameter and other stand variables. The best set of stand variables were included in the following expression, which ensures that D _m is always lower than or equal to D _g (Frazier 1981):

$$ D_{\text{m}} = D_{\text{g}} - \exp{(a_{0} + a_{1} H + a_{2} N)} $$

(9)

where D _m is the arithmetic mean diameter (cm), D _g the quadratic mean diameter (cm), H the dominant height (m), N the number of stems per hectare and a ₀- a ₂ are the model parameters. This equation was fitted by the OLS technique.

Height–diameter relationship

A height–diameter relationship was used to estimate the tree height for each diameter class. A single h–d relationship may not be adequate for all situations because it varies between stands and also with age (Curtis 1967). To solve this, the generalized h–d relationships include stand variables to localize the height predictions for each stand age.

In modelling the h–d relationship, we selected equations that are constrained to predict stand dominant height (H) from dominant diameter (D _d) and to yield 1.3 m for zero d values (to prevent negative estimates in small trees). From these, previous analyses revealed better performance of a generalized form of the model of Burkhart and Strub (1974) (Eq. 10), which was also used by Gómez-García et al. (2014b).

$$ h = 1.3 + (H - 1.3) \exp\left( (a_{0} + a_{1} H + a_{2} D_{\text{g}}) \left( \frac{1}{d} - \frac{1}{D_{\text{d}}}\right) \right) $$

(10)

where h is the tree height (m), d the tree diameter at breast height (cm), H the dominant height (m), D _d the dominant diameter (cm), D _g the quadratic mean diameter (cm) and a ₀- a ₂ are model parameters.

Practical use of the generalized h–d equation requires estimation of the dominant diameter, which was estimated from the diameter distribution, as it is difficult to project in time (Lappi 1997).

Stem taper function

When the height of the average tree in each diameter class is known, the merchantable volume from the stump to a fixed top diameter limit can be estimated directly by using a tree volume ratio equation or by integration of a stem taper function (Castedo Dorado et al. 2005; Diéguez-Aranda et al. 2006). According to these authors, the latter is the most commonly used approach.

In this study, we used the stem taper function fitted in Arias-Rodil et al. (2015b), which is based on the Kozak (2004) model. Assuming that no additional measurements will be available to calibrate the model, we considered the parameter estimates corresponding to the fixed-effects model, which was fitted by OLS:

$$ d_{\text{i}} = a_{0} d^{a_{1}} h^{a_{2}} x^{(b_{1} (h_{\text{i}} / h)^ 4 + b_{2} (1/\exp(d/h)) + b_{3} x^{0.1}+ b_{4}(1/d)+ b_{5} h^{w} + b_{6} x)} $$

(11)

where d _i (cm) is the top diameter at stem height h _i (m), d the diameter at breast height (cm), h the total tree height (m), w = 1−(h _i/h)^1/3, $x = \frac {w}{1 - (1.3/h)^{1/3}}$, and a ₀- a ₂ and b ₁- b ₆ are the model parameters.

Equation 11 cannot be analytically solved for h _i and cannot be directly integrated to obtain the volume to a top diameter limit (v _i), which implies that numerical procedures should be used. The h _i at which a specific diameter (d _i) is reached was obtained by the optimize function of R (R Core Team 2015), which uses the method of Brent (1973), considered appropriate for one-dimensional optimization. On the other hand, the total and merchantable tree volumes (v and v _i, respectively) of the average tree of each diameter class were estimated by numerical integration with the integrate function of the same software. The products of these tree volumes and the corresponding number of stems were finally aggregated to estimate total (V) and merchantable (V _i) stand volumes.

2.3 Overall evaluation of the model

Each submodel of the dynamic model was evaluated graphically and numerically. We visually inspected plots of residuals against estimated values, to evaluate the presence of heteroscedasticity, and graphs of the predictions of transition functions Eqs. 1, 2 and 3 overlaid on the observed trajectories, to assess the adequacy of each submodel. Numerical analyses were based on the root mean square error (RMSE) and the adjusted coefficient of determination for nonlinear models ($R^{2}_{\text {adj}}$). Both of these take into account the number of parameters, thus penalizing inflated model parameterization. These error statistics are expressed as follows:

$$ \text{RMSE} = \sqrt{\frac{\displaystyle {\sum}_{i = 1}^{n}(y_{i} - \hat{y}_{i})^{2}}{n - p}} $$

(12)

$$ R^{2}_{\text{adj}} = 1 - \frac{\displaystyle {\sum}_{i = 1}^{n}(y_{i} - \hat{y}_{i})^{2}}{\displaystyle {\sum}_{i = 1}^{n}(y_{i} - \bar{y})^{2}} \frac{n - 1}{n - p} $$

(13)

where y _i and $\hat {y}_{i}$ represent the observed and estimated values, $\bar {y}$ is the average of y _i values, n is the total number of observations, and p is the number of model parameters.

Although each submodel would perform adequately, this does not guarantee that the overall dynamic growth model provides reliable results. Thus, an overall evaluation should be made. Because of the lack of an independent data set, we projected the information (H, N and G) of the first inventory (35 plots) to the ages of the second and third inventories by using the transition functions. At these ages, the projections were used to estimate the total and merchantable stand volumes, as they are closely related to the economic value of a stand. For this purpose, we applied both the stand volume ratio approach and the disaggregation system.

The total stand volume estimates and state variable projections were also evaluated in terms of the critical error (Reynolds 1984), which was computed by rearranging the statistic of Freese (1960) (14), and are expressed as a percentage of the observed mean. A critical error of 10 to 20 % is generally considered realistic and reasonable in forest growth modelling (Huang et al. 2003).

$$ E_{\text{crit}} = \frac{\tau \sqrt{\frac{1}{\chi^ 2_{\text{crit}}} \displaystyle {\sum}_{i = 1}^{n} (y_{i} - \hat{y}_{i})^ 2}}{\bar{y}} $$

(14)

where τ is a standard normal deviate at the specified probability level (τ = 1.96 for α = 0.05), χcrit2 is the value of the χ ² distribution obtained for α = 0.05 and n degrees of freedom, and the remaining variables are as previously defined.

Within the disaggregation system, we applied the Kolmogorov–Smirnov test (KS): (i) to all plot-inventory combinations, to assess the suitability of the two-parameter Weibull function for describing the diameter distribution and (ii) to all projected stands (projected as explained in the previous paragraph), to evaluate the performance of the parameter recovery approach when projecting state variables. The KS test compares the estimated and the real diameter distribution. Because the estimated distribution parameters are determined from the data, Lilliefors (1967) stated that the KS-statistic existing distribution is no longer valid and should be obtained by Monte Carlo simulation. For each plot-inventory combination or projected stand, we generated 10,000 independent identically distributed pseudo-random samples under the null hypothesis (i.e. with recovered parameters), and we computed the KS statistic for each sample. This enables approximation of the KS-statistic distribution under the null hypothesis, which was subsequently used: if the KS statistic value obtained from the comparison between the estimated and real distribution of a plot exceeds the critical value at a specified significance level (obtained from the approximated distribution of KS statistic), the hypothesis that the observations belong to a Weibull distribution of the specified parameters should be rejected.

2.4 Simulator in R

The dynamic growth model of Asturias can be included in a stand growth simulator, which enables simulation of different management schedules (depending on the timing, intensity and type of cutting). Simulation of thinning operations varies depending on whether the disaggregation system (information about diameter classes) or the stand volume ratio function is used. For implementing this simulator, the R statistical software (R Core Team 2015) is a good means of transfer for the following reasons: (i) it is commonly used by statistical analysts, (ii) it enables a testing workspace and (iii) it is easy to learn to use. Reasons (i) and (ii) ease the adaptation and expansion of the code by the research community, while reason (iii) facilitates the use of the dynamic growth model by forest practitioners. More details about thinning simulation and code structure are shown in the Appendix.

2.5 Comparison with other dynamic models

Once the dynamic growth model was developed for maritime pine in Asturias, it was compared with models previously developed for the same species in Galicia (included in Diéguez-Aranda et al.2009) and northern Portugal (ModisPinaster, Fonseca 2004; Fonseca et al. 2012). These models are based on the same state variables as the present model and include a disaggregation system (they do not have a stand volume ratio equation). The Galician model includes dummy variables in the dominant height and stand basal area transition functions, as well as in the h–d relationship, to differentiate between two ecoregions (Alía Miranda et al. 2009): coast and interior. No mortality was observed in Pinus pinaster stands in Galicia; thus, the model does not include a stem density reduction equation. ModisPinaster uses a disaggregation system based on the Johnson’s S _B distribution, a h–d relationship, and a tree volume ratio function to estimate total and merchantable volumes. We used a newly developed h–d relationship proposed by Gómez-García et al. (2015) for this species and region, instead of that considered by Fonseca (2004).

The model comparison between regions was carried out on the basis of (i) projection of state variables (H, N and G), (ii) prediction of diameter distribution and total stand volume from projected state variables and (iii) the age at which the mean annual increment (MAI) of total stand volume is maximal (optimal biological rotation age). For comparisons (i) and (ii), we applied each region-specific dynamic growth model to the first-inventory plots of the data set used in the present study, projecting their state variables to the ages of second and third inventories, and then predicting diameter distribution and estimating the total stand volume, as done for the overall evaluation of the dynamic growth model for Asturias. We subsequently calculated the RMSE and critical errors in projecting the state variables and in estimating total stand volume. For the diameter distribution, we computed the mean and variance of plot-level Kolmogorov–Smirnov statistics, obtained from the comparison between predicted (by each region-specific model) and real diameter distribution.

For comparison (iii), we combined four site indices (7, 11, 15 and 19 m) and four stem densities at 20 years (500, 900, 1300 and 1700 stems ha⁻¹), to obtain 16 example stands. To use the same initial stands for comparison, initial basal area was computed from these variables and from t equal to 20 years, using Eq. 4 (with parameters to be obtained in the fitting step). Region-specific dynamic growth models were then used to simulate the stand evolution following a no-thinning schedule. Finally, for each stand and model, the age of maximum total stand volume MAI was obtained, i.e. the time at which the biological productivity is highest or the optimal biological rotation age (Avery and Burkhart 2002, pp. 353–355).

3 Results

3.1 Dynamic model

Table 2 shows the parameter estimates obtained during the fitting step for all the submodels included in the proposed dynamic growth model for maritime pine in Asturias. Note that both height growth and stem taper functions Eqs. 1 and 11, respectively, have already been developed (no error statistics shown in Table 2). For transition functions, height and stand basal area growth functions are shown in Fig. 1 (top left and top right, respectively), while mortality curves are shown in Fig. 2 (solid grey and black lines). The predicted trajectories of these transition functions were plotted for different initial conditions (values of H, N and G) at age 20 years, overlaid on the observations.

Table 2 Parameter estimates and error statistics of the dynamic growth model developed for maritime pine in Asturias

Fig. 2

Mortality curves for maritime pine in Asturias (solid black lines) and Portugal (dotted lines), overlaid on the observed trajectories for Asturias, for different stem densities at 20 years: 500, 900, 1300 and 1700 stems ha⁻¹. The Portuguese mortality submodel depends on dominant height at that age as well as stem density at the initial age; the graph was obtained from the height growth curve that passes through 11 m at 20 years

3.2 Overall evaluation of the model

Table 3 shows the RMSE and critical errors obtained when using the developed dynamic growth model for projection of state variables and estimation of total stand volume. The latter was computed using both the stand volume ratio function and the disaggregation system. The results of applying other region-specific dynamic models to our data are also included but will be discussed later. Additionally, for comparison between volume estimation alternatives, Fig. 3 shows the RMSE values in merchantable volume estimation for 0- to 30-cm top diameter limits, obtained by using both approaches. The stand volume ratio function performed slightly better and was therefore used in comparison with other dynamic models.

Table 3 RMSE (between brackets critical error in percentage, E _crit) obtained in state variables projection and total stand volume prediction of projected stands (first five rows), and mean (variance between brackets) of Kolmogorov–Smirnov statistics of predicted diameter distributions of projected stands (sixth row), when applying the region-specific (Asturias, Galicia—coast and interior—and northern Portugal) dynamic growth models for Pinus pinaster to our data

Fig. 3

RMSE in merchantable volume (V _i, to a certain top diameter limit d _i) against top diameter limits, obtained with the stand volume ratio function and the disaggregation system

The Kolmogorov–Smirnov test showed that the Weibull function successfully explained (at a 5 % significance level) the diameter distribution of Pinus pinaster in 94.1 % of the plot-inventory combinations, and when evaluating the diameter distributions of projected stands, the moments method was accurate for approximately 80 % of the stands.

3.3 Comparison with other dynamic models

3.3.1 Transition functions

According to the results shown in Table 3, the models presented small differences in transition function predictions. The models for maritime pine in Galicia yielded the worst results for height growth prediction (0.7650 and 0.7386 m, for the coastal and interior region, respectively). They also produced the poorest results for prediction of stem density decrease. For stand basal area growth, the dynamic model developed for maritime pine in Portugal provided the least accurate predictions (2.475 m² ha⁻¹).

Some existing functions for other regions even outperformed those developed for Asturias. Thus, ModisPinaster (the Portuguese model) performed slightly better for mortality prediction than the Asturias model. Additionally, the stand basal area growth model developed for the interior region of Galicia provided slightly more accurate results than the submodel developed for Asturias.

Figure 1 (left) shows the height growth curves (for site indices of 7, 11, 15 and 19 m) predicted by the corresponding submodels of all the region-specific models (region-specific prediction corresponds to prediction of Asturian submodel in top graph), overlaid on the trajectories observed in Asturias and Galicia (because we had access to the data sets used for developing the Galician models). The age range of the data from Galicia (especially for the coastal region) is narrower than for Asturias. However, the models developed for maritime pine in Asturias and Galicia performed similarly up to age 30–35 years, although the latter model provided the lowest growth rates at old ages (also compared to the model developed for Portugal). The height growth function for Portugal yielded the smallest differences between site qualities at early ages, but the largest differences at old ages.

Figure 2 shows the observed mortality in Asturias and the curves (for values of stem density of 500, 900, 1300 and 1700 stems ha⁻¹, at 20 years) predicted by the mortality functions of the region-specific models (except those for Galicia, which have not been developed). The mortality model developed for Portugal predicted higher mortality than the corresponding function of Asturias, but the difference was negligible for low-density stands.

Figure 1 (right) shows the stand basal area growth curves (for basal area values of 10, 25, 40 and 55 m² ha⁻¹, at 20 years) predicted with the stand basal area growth functions of the region-specific models (region-specific prediction corresponds to prediction of Asturian submodel in top graph), overlaid on the trajectories observed in Asturias and Galicia (as for height growth, we had access to the data sets for Galicia). All the models showed similar trends except that developed for the coastal region of Galicia, which predicted the lowest growth rates (particularly for intermediate-old ages).

3.3.2 Prediction of diameter distribution and total stand volume

Table 3 shows the results of prediction of diameter distribution and stand volume for projected stands, obtained with each region-specific model, when applied to the Asturian stands. As observed, diameter distribution was similarly predicted by models developed for Asturias and Galicia, while ModisPinaster yielded less accurate predictions. Regarding total stand volume, ModisPinaster provided the worst estimations, with a considerable difference in accuracy relative to that yielded by the models developed for Asturias and Galicia. Of the Galician models, the model for the coastal region performed better than the model for the interior region.

3.3.3 Optimal biological rotation age

Figure 4 shows the optimal biological rotation ages obtained with each region-specific model (different lines) against the stem density of the example stand considered and for different site indices (different panels). The biological rotation age decreased as the site index increased and the stem density decreased. In the comparison between region-specific models, the Asturian dynamic model generally yielded the shortest biological rotation ages, while that of the interior region of Galicia provided the longest. We also observed that ModisPinaster yielded the highest MAI at these ages (average MAI accross S and N of 18.7 m³ ha⁻¹ year⁻¹), while the model for the coastal region of Galicia provided the lowest MAI (13.3 m³ ha⁻¹ year⁻¹). The models developed for Asturias and the interior region of Galicia performed similarly (15.1 and 15.6 m³ ha⁻¹ year⁻¹, respectively).

Fig. 4

Optimal biological rotation ages against stem density of the example stand considered (at 20 years) by site indices (different panels) and region-specific models used to obtain them (different lines)

4 Discussion

4.1 Dynamic model

Correlation between residuals of stand basal area growth and mortality functions was considered in the present study by fitting both equations simultanously. This is consistent with the findings of Gómez-García et al. (2014a), who used the same approach. In the development of the mortality function, some of the models tested considered the site index as an explanatory variable, but these were not as accurate as the model finally selected. The results of other studies have related the increase in site index to greater (Eid and Tuhus 2001; Álvarez González et al. 2004; Diéguez-Aranda et al. 2005) or lower mortality (Woollons 1998), which demonstrates that the effects on mortality are not clear.

When fitting the stand volume ratio function, several variables, such as number of stems, age and site index, were evaluated for inclusion in the allometric expression of the total volume. However, as these variables provided non-significant parameter estimates or only slight improvement, and following the principle of parsimony in model development, only stand basal area and dominant height were finally included in the fitted submodel. These variables also proved useful for explaining total stand volume in other studies (e.g. Diéguez-Aranda et al. 2009, pp. 135–137; Tewari et al.2014).

Besides being the most accurate method of estimating merchantable volume (Table 3), the stand volume ratio function is easier to apply and more efficient (from a computational point of view) than the disaggregation system because it does not require iterative procedures. However, we recommend use of the disaggregation system when specified log lengths are required by the market.

4.2 Overall evaluation of the model

The critical error values obtained for the state variables and volume are within generally accepted limits (Huang et al. 2003). The error corresponding to volume was similar to those obtained in dynamic models developed for other species (Diéguez-Aranda et al. 2006; Castedo-Dorado et al. 2007). Recently, Tewari et al. (2014) reported smaller critical errors than observed in the present study, but for shorter projection lengths (1-to-3 years compared with 2-to-5 years).

Regarding practical application of the dynamic growth model, the main limitation is that we did not consider the later effect of thinning and pruning before full occupation of the additional space made available for the remaining trees, as in other studies (Amateis 2000; Álvarez-González et al. 2010; García 2013). Although the second source of data corresponds to a thinning trial, the plots are located in only six sites, which we consider insufficient to allow development of a separate equation to include the thinning effect.

4.3 Comparison with other dynamic models

4.3.1 Transition functions

The stand basal area growth model for the interior region of Galicia yielded slightly better predictions than the model developed for Asturias (Table 3). This also occurred for mortality, for which the equation developed for maritime pine in Portugal was more accurate than that developed for Asturias. This may be because we fitted the stand basal area growth and stem density reduction functions simultaneously for the stands in Asturias, which might affect the performance of each separate equation in favour of the improvement of the whole system. However, the differences from submodels developed for other regions are very slight, and simultaneous fitting accounts for correlation between residuals in stem density reduction and basal area growth, which proved to be significant.

As seen in Fig. 1 (left), the height growth data from Galicia was lacking data for old ages, which implies that extrapolation beyond the age range of the data used in model fitting may be unreliable. The age range used for Asturias was much longer (up to 68 years), similar to that used for Portugal (up to 65 years, Fonseca2004, p. 8). For these regions, the predicted growth rates at old ages and for intermediate site qualities were similar. Moreover, the equation developed for Asturias provides a good representation of the observed trajectories in Galicia (see Fig. 1, left). Therefore, the height growth model for Asturias appears reliable, and we consider that stem analyses should be conducted in old stands in Galicia to assess whether the growth at these ages is similar in these two regions.

The mortality in maritime pine stands in Asturias (see Fig. 2) and Portugal (Fonseca 2004, p. 12) was very low, which is consistent with the stand mortality observed in Galicia (Álvarez González et al. 1999). Nevertheless, mortality was considered in the models developed for the former regions by the inclusion of stem reduction functions. This seems biologically more reasonable than assuming absence of mortality; therefore, this can be considered as a deficiency of the models developed for the species in Galicia, which explains the low accuracy in mortality prediction for Asturian plots (Table 3).

For basal area growth, the amount of data collected at old ages was not satisfactory for either Asturias or Galicia (see Fig. 1, right). However, the corresponding submodel for Portugal was developed on the basis of observations of a wider age range (up to 65 years, Fonseca2004, p. 9). Therefore, new basal area growth measurements should be carried out at old ages in future studies in Asturias and Galicia. Finally, the low growth rates of the submodel for the coastal region of Galicia at intermediate-old ages are probably due to the unequal data distribution for basal area and age classes, i.e. because no observations were measured from this age range or stands with large basal areas in the coastal region of Galicia.

4.3.2 Prediction of diameter distribution and total stand volume

ModisPinaster did not prove as reliable as the other models when predicting diameter distribution of Asturias stands. This may be explained by the fact that the latter are based on the two-parameter Weibull distribution, which has no upper boundary and a zero-value lower boundary, and the Portuguese model is based on the Johnson’s S _B distribution, which has both upper and lower boundaries (Fonseca et al. 2009) and is subsequently based on a narrower diameter range. In addition, it should be taken into account that a different recovery-parameter approach must be used in each case. Nevertheless, the values of central tendency (e.g. mean, median...) of predicted diameter distributions were very similar (e.g. mean value of median between 20.3 and 20.8 cm).

Regarding prediction of total stand volume of stands in Asturias, the model developed for Portugal was not very accurate, and we found that it was overestimated by the disaggregation system (mean bias for the whole model of −23.7 m³ ha⁻¹), given that the projections of transition functions are similar to those obtained with other region-specific models. A more detailed analysis showed that tree volume submodel of ModisPinaster predicted higher volumes than those developed for other regions, e.g. for a tree of d = 22.5 cm and h = 12.5 m, the Portuguese model yields a tree volume of 0.246 m³ while those of Asturias and Galicia predict 0.220 and 0.209 m³, respectively. This difference might be caused by the data set used to develop ModisPinaster, which presents an unequal data distribution relative to that used in the present study (see Fonseca2004, p. 6).

4.3.3 Optimal biological rotation age

The limitations found both in the transition functions and the total stand volume prediction affected the optimal biological rotation ages shown in Fig. 4. For example, the higher growth rates predicted by the stand basal area growth submodel of the interior region of Galicia, relative to that of the coastal region, explains why the former predicted higher MAI values. Moreover, overestimation of the disaggregation system of ModisPinaster explains why it showed the highest MAI values. Therefore, no meaningful findings can be extracted from the comparison of biological rotation age.

5 Conclusions

The state variable projections and total stand volume prediction for the whole-stand dynamic growth model developed for maritime pine in Asturias were within the generally accepted error limits. We recommend using a stand volume ratio equation to estimate the total and merchantable stand volume, as it is more accurate and efficient than the commonly used disaggregation system. To facilitate the use of the model by the research community and by forest managers, we included the dynamic growth model in a stand growth simulator implemented in R (R Core Team 2015), which is shown in the Appendix.

Considering the similarities and differences (probably caused by the available information used for model development) between region-specific models, future research should focus on collecting more data to balance the information for age and site quality classes and perhaps to develop a single dynamic model for the whole NW of the Iberian Peninsula. This would be consistent with the results of de la Mata and Zas (2010), who did not find sufficient evidence for subdividing Galicia into the two ecoregions currently considered.

Appendix: Thinning simulation

A real growth simulator for maritime pine should enable projecting the evolution of a stand under different management prescriptions, for which thinning simulation is needed. Depending on the method used to estimate total and merchantable stand volume, thinning can be simulated in different ways.

Given that the disaggregation system categorizes the number of stems in different diameter classes, thinning can be simulated to act on each diameter class separately. Once the percentage of stems that should be removed is specified, uniform thinning is easily applied by removing the same percentage of stems in each diameter class. Alder (1979) proposed a methodology to simulate thinning from below:

$$ n_{j}= N_{\text{bt}} L \left( F(d_{j})^{1/L} - F(d_{j - 1})^{1/L}\right) $$

(15)

where n _j is the remaining number of stems in diameter class j, N _bt the number of stems per hectare before thinning, L the low-thinning intensity, computed as 1−N _r/N _bt, with N _r the number of stems to remove, and F(d) the continuous distribution function of diameters.

When a stand volume ratio function is used, thinning can be simulated from stand variables. In our case, we implemented the thinning relation proposed by Álvarez González et al. (1999), which depends on the number of stems removed and the stand basal area and number of stems per hectare before thinning (16). The value varies according to the type of thinning: 0.35 to 0.60 for thinning from below, 1 for uniform thinning and > 1 for thinning from above.

$$ R_{\text{t}} = \frac{G_{\text{r}}/G_{\text{bt}}}{N_{\text{r}}/N_{\text{bt}}} $$

(16)

where G _r and G _bt are the stand basal area removed and before thinning, respectively, and N _r and N _bt the number of stems per hectare removed and before thinning, respectively.

1.1 R code

The dynamic model developed in the present study was implemented in a growth simulator in an R script (R Core Team 2015). It allows estimation of the stand volume for different projection ages, both with the stand volume ratio function and the disaggregation system, and simulation of management prescriptions.

The simulator function (SimulateGrowth) is based on four functions: (i) a function to initialize basal area (InitializeBasalArea) if it is not provided, (ii) a function to simulate a growth interval not affected by thinnings (SimulateGrowthInterval), (iii) a function to estimate the volume at any time (EstimateVolume) and (iv) a function to simulate thinning operations (SimulateThinning). InitializeBasalArea depends on the height growth function (ProjectHeight) to estimate the site index of a stand. SimulateGrowthInterval depends on the transition functions (ProjectHeight, for height growth function; ProjectNumberOfTrees, for mortality function and ProjectBasalArea, for stand basal area growth function) and a function to recover Weibull parameters (RecoverWeibullParameters), to estimate diameter distribution, if the disaggregation system is used. EstimateVolume depends on the following: (i) a function with the height–diameter relationship (EstimateHeight) and a function to estimate the volume using the stem taper function (EstimateVolumeAtDi, which uses the stem taper function, EstimateDi), when the disaggregation system is used and (ii) the stand volume ratio equation (EstimateStandVolume). Finally, SimulateThinning applies the thinning operations according to the aforementioned methods, depending on whether the disaggregation is used or not.

In the script, a management prescription is defined in an R data.frame (R Core Team 2015) with four variables: age (t, year), thinning relation (Rt), per-unit proportion of stems per hectare to remove (pNr) and whether the thinning is uniform or not (uniform).

The simulator function (SimulateGrowth) uses the following as arguments: a data frame with stand information (stands), a data frame with management prescriptions (management.prescriptions), a top diameter limit (di, cm) and whether the disaggregation system will be used or not (disaggregation). Stump height (hst, m), first diameter class (init.dc, cm) and width of diameter classes (width.dc, cm) are optional arguments for the disaggregation system, which by default are set at 0.1 m, 5 cm, and 5 cm, respectively. The stands data frame contains one stand per row, with initial age (t, years), initial dominant height (H, m), initial stem density (N, stems ha⁻¹) and initial stand basal area (G, m² ha⁻¹). If G is not provided (i.e. NA), the initialization function is used. The simulator function returns a data frame with the stand number (stand), alternative number (alternative), age (t, years), dominant height (H, m), number of stems per hectare before and after thinning (Nbt and Nat, respectively, stems ha⁻¹), stand basal area before and after thinning (Gbt and Gat, respectively, m² ha⁻¹) and volume before thinning, after thinning and removed (Vbt, Vat and Vr, respectively, m³ ha⁻¹).

Implementation of the dynamic model equations and the growth simulator is shown below. By way of example, we illustrate use of the simulator to generate the output for the evolution of two stands (we consider for this case d _i = 0 cm): (1) t = 15 years, H = 7 m, N = 900 stems ha⁻¹, G = 15 m² ha⁻¹ and (2) t = 20 years, H = 10 m, N = 1000 stems ha⁻¹, G = NA (not available in R terminology); for a given management prescription: (i) uniform thinning at 25 years, removal of 30% of the standing trees; (ii) thinning from below at 35 years, removal of 40 % of the standing trees and (iii) final harvest at 45 years.