Dynamic Height Growth Equations and Site Index-Based Biomass Models for Young Native Species Afforestations in Spain

: The expansion of forested areas through afforestation and reforestation is widely recognized as a highly effective natural solution for mitigating climate change. Accurately assessing the potential carbon uptake capacity of newly afforested areas requires modelling tools to estimate biomass stocks, including site index curves and biomass models. Given the unique conditions in terms of tree size, uniform spacing, and tree allometries observed in young afforestations compared to natural stands, specific tools are necessary. In Spain, over 800,000 ha has been afforested with native forest species since 1992, but specific modelling tools for these plantations are lacking. Using data from 370 stem analyses collected across an extensive network of plots in young afforestations, we developed dynamic height growth and site index models for the main native species (five pines and five oaks) commonly used in afforestation in Spain. We compared various nonlinear models, such as ADA (algebraic difference approach) and GADA (generalized algebraic difference approach) expansions. The developed site index models were then used to predict the total biomass stored in the afforestation. Our results underscore the necessity for specific site index models tailored to affor-estations, as well as the potential of the established site index in predicting biomass and carbon fixation capacity in these young forests.


Introduction
Taking into account the capacity of trees and forests for sequestering atmospheric carbon through photosynthesis, the restoration of forests and the expansion of forested areas via afforestation and reforestation are considered among the most efficient natural solutions for achieving climate change mitigation goals [1][2][3].This critical role has been acknowledged by numerous international initiatives, including the United Nations Decade (2021-2030) on Ecosystem Restoration, the 1 Trillion Tree Initiative of the World Economic Forum, the Bonn Challenge aimed at restoring 350 million ha by 2030, and the European Union's commitment to plant 3 billion trees by 2030 as part of the Forest Strategy.However, the determination of the potential capacity for carbon uptake by new afforestation is a topic of fierce debate [4].A high degree of uncertainty exists around the availability of large deforested areas suitable for these afforestation programs and the potential trade-offs with the conservation of natural ecosystems (shrublands or grasslands), as well as competing land uses such as agricultural production or water supply [4,5].Moreover, another source of uncertainty is that of the carbon capture potential of newly forested land [6][7][8].An accurate assessment of the carbon uptake potential in reforestations is required at various scales, ranging from smaller farm-owned or forest-unit scales to regional or national scales.Moreover, international institutions request periodic reports on accumulated carbon in different sinks, including young afforestations.While growth, yield, and biomass models are common tools for assessing stocks in adult and natural forests, specific modelling tools for young afforestation are still lacking [6,9].Most of the existing biomass equations and modelling tools for forest species were developed for trees larger than the minimum threshold used in forest inventories.Therefore, applying these equations to trees of smaller sizes from young reforestation could result in biased estimations of carbon stocks [10].In addition, a detailed quantification of temporal changes in the carbon stock due to afforestation is necessary, e.g., for LULUCF (Land Use, Land Use Change, and Forestry) reporting, carbon footprint compensation projects, or for monitoring Clean Development Mechanism projects [11].Hence, there is an important demand for dynamic models.
Since the seminal theory by [12], which postulated a direct relationship between stand height and standing volume, the site index-defined as the height attained by the stand at a given age-has traditionally served as an indicator of potential productivity and site quality [13].Site index curves have traditionally been constructed using height growth dynamic models [14,15], which depict the pattern of height growth with age throughout the lifespan of the tree or stand.Differences in the height growth pattern can be observed between natural forests and plantations, particularly in the initial stages of development [16,17].These differences can be due to the following: As a consequence of this differential pattern in growth, using a height growth dynamic model developed for natural stands in artificial afforestation may lead to a severe bias in the site index estimation [18].In addition, due to the commonly observed erratic growth pattern in the initial stages, a site index classification based on reference ages close to the rotation length may be of little use for very young stands [19].Therefore, specific height growth dynamic models and site index curves are necessary for afforestations, especially in the initial years after plantation.
In Spain, more than 1,200,000 ha were afforested during the 30-year period from 1992 to 2021 [20,21].The main part of this area (>700,000 ha) corresponds to the forestation program for former arable lands, carried out on private land within the framework of the European Union Common Agricultural Policy (CAP).The rest of the total forested area largely consists of afforestations on public land.Of the total afforested area, approximately only one third (400,000 ha) was planted with high-timber-producing exotic species, such as Eucalyptus sp., Pinus radiata.Pseudotsuga menziesii, Populus × hybrida, and P. nigra var.austriaca.In the remaining 800,000 ha, slow-growing native species were used, mainly consisting of five pine species (Pinus pinea, P. halepensis, P. pinaster, P. nigra var.salzmannii, and P. sylvestris) and five oaks (Quercus ilex, Q. suber, Q. robur, Q. faginea, and Q. pyrenaica), with a small percentage of other broadleaves such as Castanea sativa, Fagus sylvatica, Betula sp., and Olea europaea.Unlike previous afforestation programs, the promotion of mixtures in the plantations has been incentivized in recent decades [21].The primary goal of these afforestations with native species is the restoration of degraded and non-forested ecosystems and protection against erosion, with a secondary objective focusing on both timber and non-timber (cork and nut) production.Finally, in recent years, several thousand hectares have been planted within programs aimed at carbon footprint compensation, promoted by different industries and companies (see https://www.miteco.gob.es/es/cambioclimatico/temas/registro-huella.html,accessed on 2 May 2024).
In recent decades, different height growth dynamic models and site index curves have been developed for natural stands of the main native Spanish species, including conifers [22][23][24] and broadleaves [25][26][27][28].In addition, site index curves have been developed for the exotic species used in high-timber-producing plantations, such as Eucalyptus sp.[29], Pinus radiata [30], or Pseudotsuga menziesii [31].These site index curves have been commonly employed as the basis for constructing growth and yield models and tables that provide forest managers with necessary outputs, such as total volume or biomass stocks (see [32]).However, despite the vast area occupied by recent native species afforestations in Spain, specific tools for supporting the management of these stands, such as site index curves or models for assessing stocking biomass and CO2 fixation capacity, are currently unavailable.
The main aims of this study were (1) to develop a set of dynamic height growth models for the main native species commonly used in recent afforestations in Spain; (2) to construct a family of site index curves designed for these stands; and (3) to fit site qualitybased models for predicting the total biomass stock stored at different ages.Our main hypothesis is that, in the case of young afforestations, a site index defined as the height reached by the plantation at a reference age can be a reliable indicator of productivity and biomass stocking.

Network of Permanent Plots in Young Afforestations
Between 2019 and 2023, a network of 198 permanent plots was established in young afforestations planted with the most relevant native species used in Spain over recent decades.These species included five pine species (Pinus halepensis, P. nigra, P. pinaster, Pinus pinea, and P. sylvestris) and five oak species (Quercus faginea, Q. ilex, Q. robur, Q. pyrenaica, and Q. suber).A few plots in which the species used was Q. petraea were also identified as Q. robur.This network was established with the aim of extending a previously existing network of 178 plots, set up by the Spanish Ministry of Environment 2008-2010 in afforestations on former arable land.The network included plots in plantations with ages ranging from 4 to 30 years.The plantations were selected with the aim of embracing the existing range of ages, stand densities, and site qualities for the selected species.Additional criteria for selecting plantations were a success of over 80% and the absence of thinning.Within each plantation, we installed one plot in a location with average conditions.Plots were rectangular and had a variable size (depending on the initial plantation density), in order to include a minimum of 30-40 standing trees.At the establishment of the plot, the root collar diameter (cm), breast height diameter (cm, if the tree was over 1.30 m height), and total height (m) of each tree were recorded, together with the crown diameter in a subsample of five trees per plot.General attributes of the afforestation, such as the density of the plantation, spacing, date of the plantation, soil preparation technique, and ownership, together with ecological traits as slope, aspect, rockiness, and shrub coverture, were also recorded.Additionally, a subset of 67 plots, established by INIA between 2006 and 2008 in 20-45-year-old afforestations with P. pinea, P. nigra, P. pinaster, and P. sylvestris in the Castille and Leon and Andalusia regions, were included in the network, which finally comprised 443 plots (Figure 1).In these additional plots, the tree-and plot-level measurements were similar to those recorded in the network of permanent plots.
Of the total 443 plots, 249 (56.2%) are pure afforestations, while 179 (40.4%) are composed of two species and only 15 (3.4%) comprise three or more species.Given these mixtures, we have 654 records of species × plot combinations.Table 1 includes the general characteristics of the sampled plots and species.Plots refer to number of plots where species occurs, Area: mean plot area, N: current plantation density, Hm: mean height, RCD: mean root collar diameter, Wtot: total biomass stock.

Biomass Stock
The allometric biomass models for young afforestations, developed by [8] using trees collected from the network of plots, were used to compute the aboveground biomass of each tree within each of the plots of the network (see Supplementary Material, Table S1).The individual biomasses of the trees in a plot were then summed and upscaled to the hectare (Table 1).In the case of mixed plots, and given the regular spacing of the afforestations, we assumed that each species occupies an area proportional to its representativeness in number of stems per ha, and computed the total biomass upscaled to the hectare separately for each species.Henceforth, when the term plot is used in this article, we refer to a plot × species combination.

Stem Analysis
In a subsample of the plots, a total of 1184 trees of the selected species (with average conditions) were destructively sampled for biomass estimation (Table 2, see [8] for additional information on sampling methods).Among these felled trees, three per plot (one from each tertile of the height distribution) were selected for complete stem analysis.Unlike studies concerning site index modelling in adult stands, we did not only consider dominant trees, since the processes of dominance-differentiation in young afforestations is still ongoing and not yet clearly defined.Additional criteria for selecting the trees were (1) a root collar diameter > 3 cm and (2) the possibility of obtaining at least two section disks with a diameter > 2 cm.The final number of stem analyses of the species was 370, ranging from 9 trees in the case of Quercus pyrenaica to 62 in Pinus nigra, with a greater representation of Pinus sp.than Quercus sp.(Table 1), in terms of number of trees as well as the range of heights-ages sampled.
Once felled, the stem was divided into different sections, depending on the total height of the tree: -If the tree height < 2 m: a root collar section and one section every 0.5 m, up to an end section diameter of 2 cm; -If the tree height 2-6 m: a root collar section and one section every 1 m, up to a section diameter of 2 cm; -If the tree height > 6 m: a root collar section, breast height (1.30 m) section, and then one section every 2 m, up to an end section diameter of 2 cm.
A disk was extracted from each stem section and sent to a laboratory, where it was sanded.Each section disk was scanned at a resolution of over 800 ppm, and the number of annual growth rings was visually assessed, counted, and transformed to age.As the cross-section lengths did not coincide with a periodic height growth, we adjusted the height-age data from the stem analysis to account for this bias, using the correction proposed by [33].

Height-Age Dynamic Models
Ten dynamic height growth difference models (see Table 3), derived from the original base growth equations proposed by [34][35][36], were fitted to the stem analysis data.Among these ten equations, we included models derived using both the algebraic difference principle (ADA models) presented by [37] and the generalized algebraic difference principle (GADA models) proposed by [38].The evaluated GADA models have some of the desir-able properties for height-age dynamic models, such as polymorphism and multiple asymptotes, while ADA models are polymorphic with a single asymptote (models M1, M3, M4, and M7) or anamorphic with multiple asymptotes (model M2).Table 3. Evaluated dynamic models and base original models used to develop dynamic models.

Model
Base Model Expanse Dynamic Model )

Model Fit
Two approaches have been commonly used to fit difference equations, the dummy variable method [39] and the difference method [40].The first method considers the inclusion of both global and local (subject-specific, i.e., tree-level) parameters, while the difference method eliminates the need to include subject-specific parameters in the model, and therefore only global parameters need to be estimated.The disadvantage associated with the dummy variable method is the need to define an arbitrary observed age-height pair to predict the local parameters (commonly the age-height pair of the median section [41]), while in the difference method, the main disadvantage is the need to define a specific structure for the data set [42].In our study, since many of the young trees are small and only have two or three sections, we preferred to use the difference method because (1) if one of the sections is considered as the local indicator, very few section disks will remain for prediction, and (2) fitting a specific local parameter for each tree would result in problems of convergence (a large number of parameters if compared with the number of observed pairs).In fitting the difference method, following a preliminary fit, we chose the data structure VI.This structure of the data set considers all the possible growth intervals (non-descending and descending, non-overlapping and overlapping) between height-age pairs recorded for each tree [42].In this manner, for a tree with three observed height-age pairs, H1 − T1, H2 − T2, and H3-T3, we obtained the following six intervals: To deal with the inherent serial autocorrelation derived from using observations from the same tree, we added a continuous-time autoregressive expansion of order 2 to the error term of the model eij [43]: where eij is the jth ordinary residual of the ith tree (i.e., the difference between the observed and the estimated heights of the tree i at age measurement j), dn = 1 for j > n and is zero for j ≤ n, ρn is the n-order autoregressive parameter to be estimated, and tij − ti(j−n) is the time distance (years) separating the jth from the jth-n observations.After fitting models M1-M10 to the stem analysis data for each species, we compared the models for each species in terms of fitting statistics-root mean square error (RMSE), adjusted coefficient of determination (AdjR 2 ), Akaike's information criterion (AIC), delta Akaike's information criterion (AICd), and significance of the parameter estimates.The criteria for selecting the best model for each species were (1) all the parameters being significant (p-value < 0.01); (2) AICd < 3; (3) maximum AdjR 2 and minimum RMSE; and (4) visual agreement between the fitted curves and the growth trajectories obtained from the stem analysis.

Selection of Reference Age for Site Index Curve Construction and Models for Total Biomass
One of the practical applications of the above described height-age dynamic models is to determine a site index, defined in this case as the mean height of the plantation attained at a given age.One controversial point is the selection of a base age to which the site index will be referenced, from an observed pair of values of mean height and plantation age.Inversely, the site index calculated for a given plantation can be used to project the mean height at any plantation age, allowing the construction of growth and site index curves.Among the different criteria to select a reference age, the most commonly used are the minimizing of the relative error in predicting heights at other ages [44] and visual agreement between the constructed curves and the original stem analysis (e.g., [26]).
In our case, taking into account that one of the main uses of the site index in young afforestations is to determine and project the biomass, we proposed a comparison of the site index computed using different reference ages, in terms of the accuracy in predicting the total biomass stock.For this purpose, we first used the best height-age dynamic model selected for each species in the previous phase to compute three different site indices for each of the plots and species included in the network.The indices were computed using three different reference ages.For the pine species, which present wider age ranges in both the stem analysis and the network of permanent plots, we evaluated 10, 20, and 30 years as reference ages.For the Quercus sp., we compared reference ages of 10, 15, and 20 years.For each species, we then fitted the following power model: where Wtoti is the total biomass stock (kg ha −1 ) of the species in the ith plot, Agei is the plantation age (years), and Ni: the current plantation density (stems ha −1 ).SIji refers to the site index for the species in the ith plot, computed using the best model selected for each species in the previous phase, at each of the three reference ages j.Comparisons among the three different models for each species were made in terms of AICd, AdjR 2 , and absolute mean error |E|.Given the inherent heteroscedasticity observed in biomass models, weighted nonlinear regression was used to fit model [2].
All the statistical analyses were carried out using the MODEL and NLIN procedures in SAS On Demand for Academics.

Selection of the Best Height Growth Dynamic Model for Each Species
Table 4 shows the goodness-of-fit statistics for the ten models compared over the ten analysed species.Models M8 and M9, GADA expansions of the base models by Mc Dill-Amateis and Korf, respectively, only attained convergence for one or two species, while model M10, another GADA expansion of the model by Korf, though converging for nine out of the ten species, always resulted in the worst fit among the models.The rest of the models showed quite similar values of AdjR 2 (not shown) and RMSE, pointing to the importance of using additional criteria, such as AICd, and the significance of parameters for selecting the best model.Note.RMSE: root mean square error (m), sig: level of significance ***: all the parameter significant at p-value <0.05; **: one parameter non-significant (p-value> 0.05); *: two parameter non-significant (p-value> 0.05).NC: not converge.In bold, the selected model according with the defined criteria.
After comparing all the fitted models using the proposed criteria, the best models selected for all the species-except for Quercus faginea-are based on the original base model by Chapman-Richards.Model M1, an ADA expansion based on solving the original Chapman-Richards base model by parameter α3, resulting in a polymorphic expression, was the selected model for P. sylvestris, P. nigra, P. pinaster, and Q. suber.Model M2, also an ADA expansion of the model obtained after solving by asymptote parameter α1 but showing anamorphic behaviour through multiple asymptotes, was the best model for P. halepensis and Q. robur.Model M6, a GADA expansion of the Chapman-Richards model, was the best model for P. pinea, Q. ilex, and Q. pyrenaica.For the latter parameter, a1 showed a p-value > 0.05, but as no other model attained AICd < 3 with a level of significance for all the parameters < 0.05, and the best model with all the parameters significant (model M2, AICd = 4) showed worse visual agreement of the curves against the original growth trajectories (see Supplementary Material Figure S1), we decided to retain model M6 as the best.For Q. faginea, model M3, the ADA expansion of the Chapman-Richards model obtained after solving parameter α2, was the model which displayed the best goodness-of-fit criteria, although, after a graphical analysis, we detected an illogical pattern.Thus, we finally selected the second best model for the species, which was M7, the ADA model by [36].Finally, two commonly used models, the M4 (Bailey-Clutter) model and the GADA expansion of Chapman-Richards in model M5, though converging for most of the species, were not selected for any of the analyzed species.Table 5 shows the fitting parameters for the best model selected for each species.

Site Index Curves and Total Biomass Models for the Different Species
Irrespective of the proposed reference age, the site index showed a high total biomass prediction capacity when incorporated into power age models which already included current plantation density (Table 6), reaching AdjR 2 values of over 0.82 and an absolute error |E| over or below 10,000 kg ha −1 .Given the structure of model M2, no differences due to the reference age were observed for P. halepensis and Q. robur.For the rest of the species, except for P. nigra, the best results were obtained for site indexes computed using the oldest reference age (20 years for Quercus sp. and 30 years for P. sylvestris, P. pinea, and P. pinaster).For P. nigra, the best results were obtained for a reference age of 10 years.However, since the species AICd for a reference age of 30 was <1, we decided to select 30 years as the reference age for pine species and 20 years for oak species.Table 7 shows the parameters of the total biomass model (Equation ( 2)).Site index curves were drawn up using these reference ages (Figures 2 and 3).Table 6.Goodness-of-fit statistics for the models for total biomass (Equation ( 2)) as a function of the reference age selected for defining the site index.

Discussion
In this study, we have developed height growth dynamic models and site index curves for young afforestations of the ten main native forest species commonly used in plantations in Spain over recent decades.The models derived from Chapman-Richards performed best for nine out of the ten species evaluated, highlighting the flexibility of the original model in adapting to different growth patterns.Surprisingly, there was no clear preference for the GADA expanded equations over the simpler ADA models, as the former were only selected for three of the species (Q.pyrenaica, Q. ilex, and P. pinea).One possible explanation for this outcome could be that, in the case of young afforestations far from rotation age, one of the main advantages of GADA models-the ability to consider varying asymptotes in polymorphic growth [38]-may not be as relevant.However, the erratic behaviour observed in the younger stages of the trees, along with the potential effects of pre-plantation labour [19], underscore the necessity for polymorphic models.Afforestations in Spain have traditionally employed mechanical soil preparation methods such as tilling, subsoiling, ridging, and/or digging individual holes using excavators [21].Over recent decades, the plantation of container-grown seedlings hardened in nurseries has been preferred over sowing [45], and initial post-planting techniques, such as weeding or tree shelters, have commonly been used.Finally, many of the afforestations have taken place on former arable land, where previous labours and fertilization have led to differential physical and chemical properties of the soil [46,47].Consequently, differential patterns of height growth, especially at younger ages, are expected between afforestations and natural stands [17].As an example, Figure 4 compares the site index models for young afforestations of P. pinea and Q. suber, fitted in the current work, with the site index curves for natural stands of the species commonly used in Spain [22,26].As expected, plantations exhibit an accelerated height growth in the initial stages compared to natural stands, especially at the most productive sites.Applying the site index models for natural stands to afforestations, especially during the initial stages, can lead to a clear overestimation of site productivity, severe biases in biomass or volume estimation, and inaccurate stand prescriptions for plantations [19].Therefore, specific models tailored for young afforestations are essential.The selection of the reference age and validation of the constructed height growth models were conducted by assessing the suitability of the defined site indexes to predict the biomass stocking of the young afforestations across a broader network of permanent plots.Little differences were observed between the three reference ages evaluated for each species, although older reference ages (20 years for oaks, 30 years for pines) commonly resulted in more accurate biomass estimates.These site index-based biomass models allow for a reliable estimation of the actual biomass in afforestations at a given plantation age, by simply sampling the mean height and current density, which can be achieved through field-sampling inventories or the use of remote-sensing techniques such as airborne LIDAR or aerial photography [48].In this regard, the proposed models are of great utility for the ex-post (actual) estimation of biomass stocks in afforestation projects [49].Moreover, these models enable the evolution of the biomass stocks over time to be forecast, facilitating ex-ante (projected) estimates of change in biomass and carbon stocks.Both ex-ante and ex-post estimates of biomass and carbon stocks and changes are crucial for monitoring and certifying afforestation-reforestation projects within the framework of Clean Development Mechanisms [11].In the case of multi-species plantations, the site index and biomass stock should be separately computed for every species and then the total biomass per ha computed, assuming that each species occupies an area proportional to its representativeness in the number of stems per ha.
The constructed site index curves and site-based biomass models serve as valuable tools for comparing plantations with the same or different species at various ages, in terms of expected biomass stocks.In this regard, the combined use of the constructed models with environmental niche models, which define the suitability of different species for ecosystem restoration [50], provides a highly useful tool for helping managers and landowners to decide which species to include in plantations.As site index serves as a proxy or an aggregated index of different environmental factors, future research efforts should focus on identifying the edaphic, physiographic, and/or climatic factors that define site productivity in afforestations and construct site index models based on these attributes [28,51].

Conclusions
In the present work, we have used a set of 370 stem analyses to construct a set of height growth dynamic models for the ten main native species commonly used in afforestation in Spain.The dynamic models were fitted using the difference equation approach and were therefore used to develop site index curves.Models for predicting the aboveground biomass stock of young afforestations of the studied species were then constructed, including as predictors the site index, current density, and plantation age.Height growth dynamic models, site index curves, and biomass stock models are useful tools for projecting the evolution of an afforestation and its biomass stock over time, and for classifying and comparing afforestations of different ages.

Figure 1 .
Figure 1.Net of permanent plots in young afforestations in Spain.Dark dots correspond with pine afforestation, clear triangles with oak afforestations.

Figure 2 .
Figure 2. Site index curves (solid lines) for the main native pine species.The curves are forced to pass at a reference age of 30 years through heights 4, 8, 12, and 16 m for P. sylvestris, P. nigra, and P. halepensis; 6, 11, 16, and 21 m for P. pinaster; and 4, 7, 10, and 13 m for P. pinea.Dotted lines represent the growth trajectories for each of the analyzed trees.

Figure 3 .
Figure 3. Site index curves (solid lines) for the main native oak species.The curves are forced to pass at a reference age of 20 years through heights 3, 6, 9, and 12 m for Q. robur-petraea and Q. pyrenaica;

Figure 4 .
Figure 4.A comparison between the site index curves for young afforestations (solid line) for P. pinea and Q. suber constructed in this work, with the previously existing site index curves for the species in natural stands (dotted lines) by Calama et al. (2003) [22] and Sánchez-González et al. (2005) [26].

Table 1 .
Mean characteristics of sampled plots.

Table 2 .
Mean characteristics of stem analysis used to fit height growth dynamic models.

Table 5 .
Parameter estimates for the best height growth dynamic model selected for each species.

Table 7 .
Parameter estimates for the total biomass model as a function of age, plantation density, and site index, computed using, as a reference age, 30 years for pines and 20 years for oaks (Equation (2)).