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.

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.