Introduction Top

The use of growth and mortality models is essential in forest management to predict the development of forest stands under various silvicultural systems, assisting forest owners and managers in decision-making at both stand and forest levels (Bolte et al., 2009). As forest management practices continue to intensify, the need for these models becomes increasingly apparent (Crecente-Campo et al., 2010).

Growth and yield models are commonly categorized as either stand-level or individual-tree-growth system (Vanclay, 1994). An individual tree growth model is composed by different equations (Burkhart & Tomé, 2012). These equations have been classified in this work into three categories: (i) auxiliary functions to estimate missing variables; (ii) increment, growth, and yield models; and (iii) regeneration and mortality equations. In the application of growth and yield models, auxiliary functions are often required to estimate missing variables. This is because some variables, such as crown width (cw), canopy base height (cbh) or total height (ht), may only be measured on a subset of trees or because certain input variables cannot be directly measured. Increment, growth, and yield forest models are developed to quantify the growth of a forest and are commonly employed to predict the future condition of a forest and to assist in evaluating alternative silvicultural guidelines (Vanclay, 2006) . Both tree mortality and probability of natural regeneration play an important role in forest dynamics. Tree mortality can be driven by competition indices, site variables, and can ultimately lead to self-decline of the forest (Monserud & Sterba, 1999; Bravo et al., 2011). Recruitment models are used to predict the number of new trees that meet a specified minimum size threshold, as well as their biometric properties such as diameter at breast height (dbh) and height. These models are typically calibrated to ensure that the expected stand structure and canopy species composition are consistent or by utilizing empirical regeneration data (König et al., 2022).

The Aleppo pine (Pinus halepensis Mill.) is naturally distributed in the Mediterranean area of Spain (del Río et al., 2008). It covers approximately 7% of the forested land in Spain (Herranz, 2000), with 15% of the plantations established between 1940 and 1980 consisting of this species (Vadell et al., 2016). In general, plantations differ from natural stands mainly in the number and distribution of seedlings at early stages (Murray, 1986). According to Vuokila (1980), conifer plantations, when managed under the same silvicultural guidelines as natural stands on the same site, usually do not outperform in volume production. However, in practice, natural and artificial stands are hardly comparable, since natural stands typically have a higher tree density in the early stages, while plantations usually have a tree density of less than 1500 trees/ha. This difference makes it difficult to apply the same silvicultural guidelines to both, and often results in greater timber production in plantations. Though Cabanillas-Saldaña (2010) conducted a growth and yield model for natural stands of Aleppo pine in Aragón (Spain), no set of equations has been tailored for this species’ plantations in Spain.

This study aims to create a model system for simulating the development of P. halepensis plantations in Aragón (Spain) at the individual tree level. The system includes nine different models. Finally, based on the demands of forest managers in the area, and using this set of models, we have simulated different forest management scenarios (synthesized in three yield tables) represented by: timber production, soil protection and biodiversity enhancement.

Material and methods Top

Study area

The study focuses on forests of Aleppo pine sited in the Aragón region (Fig. S1 [suppl]), Northeast Spain. The pine stands studied were generally afforested between 1925 and 1990. These stands tend to present less structural complexity and poor biodiversity (Granados et al., 2016) compared to natural stands, one extra reason encouraging us to develop a differentiated model for plantations.

Field data

Two databases from permanent plots have been used in this work: (i) Second and Third Spanish National Forest Inventories (SNFI), https://www.miteco.gob.es/es/biodiversidad/temas/inventarios-nacionales/inventario-forestal-nacional/default.aspx, and (ii) plots established ad hoc by the Diputación General de Aragón (DGA. The SNFI plots were only used to construct 3 of the 9 models, as they provided field data that were too expensive to measure in the DGA plots. A subsample of 379 trees measured exclusively in 104 plots of the Second SNFI was used to characterize some canopy parameters and the diameter under bark (dub). The DGA plots were used to develop almost all the models because they covered the entire site index and stand age ranges in the area of interest. The 58 field plots from the DGA dataset were collected in April 2017, where a total of 1684 trees were identified and measured with two cross diameters. The number of height measurements per plot varied from 8 to 10, depending on the plot size, resulting in a total of 515 measurements used to model crown base height (cbh) and height-diameter ratio or slenderness coefficient (hd). Additionally, cores were extracted from six trees per plot (348 trees). The extracted cores were measured with a measuring stage (Velmex Inc., Bloomfield NY, USA) under magnification to the nearest 0.01 mm. A natural mortality model was constructed with the DGA data using a backdating procedure (more detailed information about this procedure at Hann & Hanus, 2001) to rebuild stand attributes of each plot 10 years ago. Only 16 out of the 58 plots measured showed mortality in the last 10 years. A total of 1,720 trees were used (1,684 live trees and only 36 dead trees), indicating an average mortality of 2.09% over the 10-year period. Finally, 48 trees were felled to construct a stem taper function. Table S1 shows the abbreviations used for the variables in the field data, while Tables S2-S4 show the main descriptive statistics for these datasets.

Model fitting

A total of nine models were fitted, representing the three categories: (i) auxiliary functions for estimating missing variables (Eqs. 1 to 5), (ii) increment, growth, and yield models (Eqs. 6 and 7) and (iii) regeneration and mortality equations (Eqs. 8 and 9). The accompanying Appendix on the Forest System’s website provides a comprehensive explanation of the methodology employed in constructing the models, offering additional in-depth details. Table 1 shows both the mathematical formulation of the models and the database used to construct each of them. Model 8 (occurrence of regeneration) does not have an algebraic expression, as we used decision trees to predict the output of the three strata analyzed (pure stands of Aleppo pine, pure stands of Quercus spp and mixed stands of both). A logistic model to predict the probability of regeneration was first fitted, but resulted in a trivial model where all plots had the same probability of regeneration. For that reason, we shifted to a non-parametric modelling procedure. Regarding Eq. 9, we used a stochastic procedure to apply the mortality model. Thus, the probability of survival calculated for a given tree is interpreted as the percentage of the trees per hectare with the same diameter that will survive to the next 10-year period. We also explored the deterministic alternative, but results predicted extraordinary rates of mortality in high stocking stands, which severely contradicted what field data showed.

We employed SAS® software to fit the majority of the models (SAS Institute, 1999). The least squares method was used for parameter estimation of the linear models, utilizing the REG procedure of SAS/STAT software’s capabilities. For fitting nonlinear models, the MODEL procedure of SAS/ETS software’s capabilities was used. The mortality equation was fine-tuned via the weighted maximum likelihood method, utilizing the LOGISTIC procedure of SAS/STAT software’s capabilities. Finally, the ingrowth model was constructed using decision trees with the R packages rpart (Therneau et al., 2019). All parameters in each model were checked for biological meaning; otherwise, the model was discarded. Root mean square error (RMSE) and the coefficient of determination (R²) were used as goodness-of-fit statistics. In the case of the mortality model, the value of the percentage of concordant pairs and the area under the ROC (receiver operating characteristic) curve were used as goodness-of-fit indicators. For the occurrence of regeneration model [8], the success rate of presence/absence of regeneration was used.

Simulation of different forest management scenarios

Different forest planning scenarios were evaluated by simulating several forest conditions based on a range of yield classes, initial stocking, number and intensity of thinning, rotation, and age for first thinning. For this purpose, the models developed in this study were used as shown in Fig. 1, where from initial conditions (t₀) we could estimate the final conditions after 10 years (t₁₀). All scenarios are also presented in Tables S5-S15 [suppl]. The three most representative scenarios are shown in Table 2 where Timber Production scenario (TP_4TH_30/15) is characterized by four intermediate thinnings with a rotation period of 15 years and the first thinning at age 30; Soil Protection (SP_4TH_45/10) is characterized by four low thinnings with a rotation period of 10 years and the first thinning at age 45 and; Biodiversity Enhancement (BE_3TH_45/20) is characterized by three intermediate thinnings with a rotation period of 20 years and the first thinning at age 45. In all three cases, the same site index (SI) has been considered (SI=14 meaning Ho=14 m at age 60).

Figure 1. Workflow to simulate forestry scenarios from the models developed. Variables are defined in Table S1 and models in Table 1.

Results and discussionTop

The estimated parameters were statistically significant (p < 0.05). All models explained a large part of the variability of the data (R² ranging from 0.71 to 0.99) with minimal errors. Figs. S2-S4 [suppl] show the plots of the residuals versus the predicted values of all the developed models. In all cases, the basic assumptions of the regression models were achieved. Only in the case of models 5 and 6, to avoid autocorrelation of the residuals, a CAR(1) and AR(3) error structure was used, respectively (Hirigoyen et al., 2021).

Goodness of fit obtained in this research agrees with the results obtained by other authors. As expected, bark thickness [1] increases with increasing diameter at breast height (dbh). Total height [2] increases with increasing dominant height (Ho), dominant diameter (Do) and dbh. In this case, Cabanillas-Saldaña (2010) on the same species but on natural stands in Aragon, obtained slightly worse results (R² and an RMSE of 0.83 and 0.942 respectively). Crown width [3] increases with increasing dbh but decreases with increasing ht. Height to crown base [4] gets bigger when the relative spacing is larger (expressed as the RS index), other factors being equal. Moreover, the model predicts higher hcb in stands on better yield classes. Similarly, trees under strong competition (higher balmod values) will present shorter crowns. The stem taper function [5] depends on dbh and stem slenderness (ht/dbh). Cabanillas-Saldaña (2010) obtained slightly worse results (RMSE=1.74 and R²=0.95) than us (RMSE=1.18 and R2=0.97). The best backdating growth model [6] found depended on Hart-Becking index (RS), basal area (G) and SI. We found that ten-year diameter growth fluctuated between 0.75 and 14.0 cm depending on yield class, competition, and tree size. The difference between dub₂ and dub₁ was greater the larger the relative spacing, the smaller the basal area and the better the yield class. The driving variables of the individual-tree diameter growth [7] were dbh, G, SI and balmod. When yield class increased, then diameter increment over the last 10 years also increased. The model [8] estimates the occurrence of natural regeneration through variables that are usually readily available in forest inventories. Confusion matrix shows than in total, 75.76% of the stands were correctly classified, being more accurate in predicting the absence of regeneration (83% success rate) than its presence (66% success rate). Information related to the three models [8] developed is shown in Fig. S5. We obtained G, balmod, and RS as driving variables for mortality. If G or balmod increases, the probability of survival decreases. Conversely, when the RS increases, the competition is reduced and consequently the probability of survival increases.

Figs S3-S5 [suppl] present a graphical inspection of the magnitude and distribution of residuals for all possible combination of variables. Fig. S6 shows a graphical comparison between the density management diagram created by Cabanillas-Saldaña (2010) and the “TP_4TH_30/15” scenario. We can observe the different pattern of volume and Dg evolution over time in both cases. In natural stands, the evolution is smoother from the beginning to the end of the rotation (Dg_initial =12.7 and Dg_rotation =27.4cm; V_initial = 70.2 and V_rotation=177.9 m³·ha^-1), while in the case of plantations, these values range between 11.3 and 33.1 cm for quadratic mean diameter (Dg), and between 40.7 and 275 m³·ha^-1 for volume.

In order to facilitate a user-friendly environment for their use, all the models of the system were integrated into a web-based application SIMANFOR, under the name of “Phragon_2017” (acronym for P. halepensis reforested in Aragon).

The implementation of this tailored individual tree growth system, designed specifically for even-aged Aleppo pine plantations in Aragon, significantly enhances the efficiency and effectiveness of forest management practices in the region. This tool holds great promise for assisting forest managers in their future endeavors.