Species Mixing Effects on Height – Diameter and Basal Area Increment Models for Scots Pine and Maritime Pine

Models that incorporate known species-mixing effects on tree growth are essential tools to properly design silvicultural guidelines for mixed-species stands. Here, we developed generalized height–diameter (h-d) and basal area growth models for mixed stands of two main forest species in Spain: Scots pine (Pinus sylvestris L.) and Maritime pine (Pinus pinaster Ait.). Mixed-effects models were fitted from plot measurement and tree rings data from 726 Scots pine and 693 Maritime pine trees from mixed and pure stands in the Northern Iberian Range in Spain, with the primary objective of representing interactions between the species where they are interspersed in mixtures of varying proportions. An independent dataset was used to test the performance of the h-d models against models previously fitted for monospecific stands of both species. Basal area increment models were evaluated using a 10-fold block cross-validation procedure. We found that species mixing had contrasting effects on the species in both models. In h-d models, the species-mixing proportion determined the effect of species interactions. Basal area growth models showed that interspecific competition was influential only for Maritime pine; however, these effects differed depending on the mode of competition. For Scots pine, tree growth was not restricted by interspecies competition. The combination of mixed-effect models and the inclusion of parameters expressing species-mixing enhanced estimates of tree height and basal area growth compared with the available models previously developed for pure stands. Although the species-mixing effects were successfully represented in the fitted models, additional model components for accurately simulating the stand dynamics of mixtures with Scots pine and Maritime pine and other species mixtures require similar model refinements. Upon the completion of analyses required for these model refinements, the degree of improvement in simulating growth in species mixtures, including the effects of different management options, can be evaluated.


Introduction
Despite evidence that many mixed-species forests offer greater potential to supply ecological and socio-economical goods and services than monospecific forests [1], quantitative silvicultural study area has a mean annual temperature of 9.0 • C and a mean annual precipitation of approximately 801 mm.Plots are located at an elevation ranging from 1090 to 1277 m a.s.l.Soils are acidic (pH 3.9-5.4)with sandy loam to sandy texture, low cation exchange capacity (2.4-18.1 cmol c kg −1 ), and medium to low water-retention capacity (1.5-18 g cm −2 ) [35].Each triplet consisted of three circular plots of 15-m radius within one km of each other, including two plots dominated either by Scots pine or Maritime pine and one mixed plot that contained both species with a mixing proportion of at least 75%-25% of the total basal area.Plots within triplets had similar site conditions, age, and density, and belonged to the same management compartments with the same silviculture regime (see details in [33]).The stands were at approximately carrying capacity, (60% to 100+% of known maximum density), and none of the plots had been thinned for at least 10 years.Both species were in the same age phase comparing pure and mixed-species stands by triplet, ranging from mature (45 to 50 years) to old stands (120 to 140 years).Site quality at age 100 years indicated moderate to low growth conditions according to the top height curves of the pure stands developed for Scots pine [29,30] and Maritime pine [28].Some variations in stand age and site conditions were tolerated among triplets to cover the stand variability in the study area.
All stems >7.5 cm in diameter were inventoried and measured for diameter at breast height over bark (d, nearest 0.1 cm).For each tree, total height (h) and height to crown base (hcb, base of the lowest live primary branch) were measured using a vertex hypsometer to the nearest 0.01 m.Stand variables were calculated from measurements of triplets (Table 1), including dominant diameter (Do, cm) and dominant height (Ho, m) based on the 100 thickest trees per hectare of the stand, quadratic mean diameter (dq, cm), and stand basal area (BA, m 2 ha −1 ).d: diameter at breast height (cm); h: total tree height (m); BAI: basal area periodic annual increment (cm 2 5 years −1 ); CR: crown ratio; Rh: ratio of tree height and stand dominant height; Ho: stand dominant height, average h of largest 100 trees per ha by d (m); Do: dominant diameter of 100 largest trees per ha by d(cm); dq: quadratic mean diameter (cm); d/dq: tree social status index.SI: site index as the dominant height (m) at age 100 for Scots pine [29,30] and Maritime pine [28]; m: mixing proportion of target species.BA: stand basal area (m 2 ha −1 ); BAL: basal area of larger trees (m 2 ha −1 ); subscripts inter and intra designate interspecific and intraspecific BA or BAL.
Increment cores that included at least the last 15 years were extracted at a height of 1.30 m, avoiding dead or very suppressed trees.A total of 328 and 408 Scots pine from mixed and pure stands, respectively, and 234 and 459 Maritime pine trees from mixed and pure stands, respectively, were sampled.All of the cores were mounted, sanded until tree-rings boundaries were clearly visible, and scanned at 1800 ppi image resolution.For each cored tree, tree ring widths (mm year −1 ) were dated and measured from the scanned images using the measuRing R-package [36].We performed cross-dating using species-specific marker years, narrow rings, followed by statistical confirmation of quality tree-ring series synchronization using the dplR R-package [37].To fit basal area growth models, past five-year growth (2009-2014) needed to be related to initial conditions, which was referred to as 'backdating'.For all of the cored trees, diameter over bark at breast height in 2009 was computed using the diameter measured in 2014, and then subtracting the tree ring width measurements, and applying species-specific bark factors [27].Except for crown ratio (CR, ratio of crown length and tree height), all of the tree and stand variables were backdated.The ratio of tree height to stand dominant height (Rh) was backdated using the h-d function fitted in this study.

Generalized Height-Diameter Functions
In this study, we considered five functions tested previously to fit the h-d relationships for both species in monospecific stands [27].Additionally, we tested four asymptotic nonlinear functions that have been widely used for modeling h-d relationships (Table S1).All of these models were modified to include various tree-level and stand-level variables, and the parameter representing the asymptote of the models was expanded as a linear function of the selected covariates.Including stand-level variables in the formulation of h-d models helped consider the local environment in the tree height estimates and ensured that the curves depicted within-stand trends in a single year of measurement.Height-diameter relationships are influenced by the tree social status and vigor, site quality, stand developmental stage, stand density, competition, and stand composition [9,[38][39][40][41].We evaluated the contribution of stand and tree variables (Table 1) for potential improvement of the h-d functions tested, such as dominant height and plot mean quadratic diameter to reflect the stand development stage, the ratio between the diameter of the subject tree and the plot mean quadratic diameter as a social status index, and plot basal area and the basal area of trees larger than a subject tree (BAL, m 2 ha −1 ) as a competition index.Explanatory variables were selected for each model if they were significant and improved the model likelihood, as indicated by the Akaike information criterion (AIC) [42].
Given that the h-d models developed here are intended for mixed stands of Maritime pine and Scots pine, we used species mixing proportion (m) to identify the significant effects of species interactions on h-d relationships.Mixing proportions for both species were calculated based on the stand density index (SDI) and weighted by equivalence coefficients that compare species-specific growing space requirements [16,33].
Since h-d relationships are expected to vary by species [27], models were fitted by species using nonlinear mixed-effects models (NLME) to consider the hierarchical levels of the data (multiple trees in a plot, plots nested in triplets).This nested structure could have caused high correlation among observations within a sampling level, so NLME alleviated this lack of independence of error terms [42].Accounting for within-triplet correlations did not significantly improve the models based on the likelihood ratio test.A general NLME [42] can be defined as follows: where h ij is the observed height of tree j within plot i; f is a nonlinear function of the 1 × p covariate; ϕ i is a parameter vector for plot i; and ε ij is the tree error term that is assumed to be independent and normally distributed with a mean of zero and constant variance σ 2 .In matrix form: where h i is the n i × 1 observation vector for the i th plot; X i is the n i × p design matrix for the i th plot; ϕ i is the parameter vector for the i th plot; and ε i is the n i × 1 vector of residuals for the i th plot.The parameter vector ϕ i can be split between fixed and random components [α', β i ']' to form a nonlinear mixed-effects model as follows [42]: where α i is the q × 1 vector of the fixed parameters common to all trees and plots (q is the number of fixed parameters in the model), β i is the r × 1 vector β of random parameters associated with the i th plot (r i is the number of random plot parameters in the model), and A i and B i are matrices of fixed and random effect covariates, respectively, for the i th plot.We tested different combinations of random parameters to determine which parameters should include both fixed and random components.First, all parameters on fixed effect variables were assumed to be random, with a general positive definite variance-covariance structure for the random effects.If this model failed to converge, then the number of random parameters was reduced to achieve convergence.Non-nested fitted models were compared by AIC after fitting by maximum likelihood (ML).In addition, we screened for large correlations between estimated random-effects parameters to evaluate whether the NLME model was over-parameterized.The best-fitting generalized h-d functions were re-estimated using the unbiased restricted maximum likelihood method (REML) with R-package 'nlme' [42].From the total generalized h-d functions fitted with NLME, Equation (4) was concluded to best describe the variations in the h-d relationship for Scots pine and Maritime pine, based on the lower AIC and RMSE values (Table S2): where h ij is the height of tree j in plot i, β 0 to β 3 and α 1 are considered fixed parameters; u i is a random plot effect parameter; ε ij is an error term that is assumed to be independent and normally distributed; and m is the species mixing proportions of target species, m PS for Scots pine and m PT for Maritime pine models; all of the other variables were previously specified in Table 1.Convergence was possible when only one or two parameters of the model were treated as random.However, most of the models with two parameters specified as random showed high correlations between the random parameter, indicating over-parameterized models that can produce poor estimates of the standard errors on the fixed-effect parameters in the model [42].For both species, the model treating α 1 as a random parameter at the plot level gave the smallest AIC value.Graphical analysis showed within-sample plot heteroskedasticity in the residuals, even after the inclusion of the random effects.The power variance function in Equation ( 5) was included in NLME h-d models; it homogenized the residual variance most effectively among three-variance functions evaluated, including exponential, power, and constant-plus-power functions [42]: where ε ij are the residuals at the lowest level of grouping (plot); σ 2 ε is the initial variance for the residual; δ is a parameter to be estimated from the data; and x ij is the observed tree height.

Tree Growth Model
We fitted the basal area increment model using the same NLME approach, which was based on a potential growth function multiplied by a modifier function [43,44].The potential growth estimates the maximum growth depending on tree diameter, which is then reduced by an amount proportional to competition or other limiting factors [43].The function is formulated as: where BAI is the predicted five-year basal area increment of a tree (cm 2 5 years −1 ); d is the initial diameter (cm); α 1 and α 2 are parameters in the potential growth portion of the equation to be estimated from the data; β 0 , . . ., β k are the parameters in the exponentiated linear modifier function with k explanatory variables (V 1 , . . ., V k ); and ε is the residual error.We considered CR as a covariate to indicate tree vigor; Rh to indicate competitive status of tree in the stand; and SI, or site index, as a measure of stand productivity (Ho expected at 100 years breast height as estimated from species-specific site index curves [28][29][30]).We assumed that SI in the pure stands of the triplet was a valid estimate of the site index of both species in the mixed-species plot.Finally, we included a competition term (COMP) as a combination of both a one-sided (asymmetric) and a two-sided (symmetric) competition used to adjust the potential to actual growth.We compared BA versus SDI to determine the most suitable explanatory variable for competition.Both metrics were distance-independent measures of competition that do not require tree spatial coordinates, because tree locations on the plot are not available in most of the stand inventories in the region.First, we used the basal area of trees larger than the subject tree (BAL) as the size-asymmetric competition, and the total stand basal area (BA) as the size-symmetric competition.Further, SDI was expressed in relative terms (SDIR), as the ratio between the measured SDI and the maximum potential stand density index for each species (SDI max ) [45].We used SDIR to express the size-symmetric competition and relative stand density index of trees larger than the target tree (SDIRL) [11,46] as an indicator of size-asymmetric competition for light.Parameters to estimate SDI max depend on the species, and were taken from species-specific maximum stand density functions [34,47].Covariates were included in the following order: crown ratio, competition structure, Rh, and site index.Only significant variables were included, and we used AIC and a likelihood ratio test to evaluate the overall contribution of the added covariate relative to a simpler model.
To evaluate the influence of species-mixing on the basal area growth model, competition indices for each species were calculated by pooling both species together and then splitting them into intraspecific and interspecific components [46].Thus, the resulting competition structures might express the adverse, neutral, or even positive interspecific competitive effects, e.g., when parameter estimates for interspecific competition are higher, non-significant, or lower, respectively, than the intraspecific variable.Then, these new variables were used to test whether distinguishing between intraspecific and interspecific competition results in a better fit of the basal area growth than ignoring the species contributing to the total symmetric and asymmetric competition.We identified the best competition structure taking into account model parsimony using an information-theoretic approach [48].The model with the lowest AICc (Second-order Akaike Information Criterion) and greater Akaike weight (w i ) was considered the best and most parsimonious model [48] for the observed data relative to the set of alternative models with different competition structures (see Appendix S2 in Supplementary Material).
Given the hierarchical nature of the data, we used NLME to fit a model for predicting the basal area growth while accounting for possible autocorrelation among trees within plots and plots within triplets.The significance of the random plot and triplet effect was based on the likelihood ratio test (p < 0.05) between nested models in the presence of all non-collinear fixed predictors [42].For both species, the random triplet effect was non-significant.We first included plot as a random effect in the α 1 parameter of the potential diameter growth component, and then in the intercept (β 0 ) of the modifier function, but the latter step did not improve the model performance (the intercept was not significant).Equation ( 7) depicts the base nonlinear mixed-effects model that was used for predicting the tree basal area increment of Scots pine and Maritime pine in mixed and pure stands: where BAI ij is the basal area increment of tree j within plot i; α 1 , α 2 , and β 1 are parameters to be estimated from the data; COMP ijk represents the k th competition term experienced by tree j in plot i according to the different structures described above; β k are parameters for competition effects estimated from the data; u i is the random plot effect; and ε ij is the error term.Random effect u i and error term ε ij are assumed to be normally distributed with a mean of zero and a variance σ i 2 and σ ε 2 , respectively.Fitted models showed non-constant variance in the residuals for both species.This problem was addressed by weighting observations using the same variance function as in Equation ( 5), testing d, BAL, and BA to model the residual variance.We verified mixed-model assumptions graphically (quantile-quantile and residual plots).Individual-tree basal area increment models were fitted by species using the nlme R-package [49].ML was used to compare different competition structures, and then REML was used to fit the parameter estimates and variance components for the final model selected for each tree species.

Model Evaluation
Predictions of the h-d function and basal area increment model for both species were examined to assess their predictive ability.We used the Third Spanish National Forest Inventory data (NFI) as an independent data source to quantify the model errors of the h-d generalized functions (Figure S3).A total of 213 permanent plots of mixed and pure stands from the NFI network located in the Northern Iberian Range were selected for this purpose, including 47 mixed plots, 109 pure plots of Scots pine, and 57 pure plots of Maritime pine.Plots were considered pure stands when the percentage of the total basal area of the target species was higher than 90%, and mixed stands when both species together accounted for at least 90% of the total basal area, and the proportion of each target species in the mixture plots was higher than 15%.Predictions from the most parsimonious generalized h-d functions that expressed species-mixing effects in their structure were compared with the predictions of the models proposed for the same species parameterized with data from monospecific stands [25,26].
There were some limitations with using the NFI data to evaluate the basal area increment models.Stand age and height to crown base were missing, so the CR could not be calculated.Instead, we evaluated the predictive quality of the basal area increment models using a 10-fold block cross-validation procedure.The original dataset was split at the plot level into 10 subsamples of approximately equal size.The final models were fitted 10 times for each plot, with one of the subsamples omitted at each fit for training the model.Each of the 10 fitted models was tested with the observation of the omitted subsample, and statistics were computed from the differences between the predicted (P i ) and observed (O i ) values across all plots.These predictions provided a pseudo-independent measure of the mean bias, root mean square error (RMSE), and model efficiency (EF), which were computed as follows: where O i and P i were the observed and predicted values, respectively, for the ith observation; n was the total number of observations; and O was the mean of all the observations of the variable.In Equations ( 4) and ( 7), each predicted value P i was estimated by setting µ i set to zero; i.e., predictions were based only on fixed effects from the final models; any bias in the model that was attributable to the expectation of a model with nonlinear random effects not including the population mean condition fixed effects was assumed to be minor [38].

Generalized Height-Diameter Selected Models
After the pre-selection of competing variables and comparison among candidate models (see models ranking based on AIC and RMSE values in Table S2), Equation ( 4) was the most plausible model to describe h-d relationships in Scots pine and Maritime pine.Dominant height and quadratic mean diameter were the most important variables in the final models; further total basal area and tree social status index were non-significant, and did not improve model performance.The inclusion of the mixing proportion of the target species improved the fit statistics of selected equations, reducing the AIC value for Scots pine by 7.2 and for Maritime pine by 1.7 (Table 2).The height predictions of the models covered most of the observations and showed adequate prediction accuracy (Figure S4).On inspecting the residuals versus predicted height from the best models for each species, we observed a homogeneous variation of the residuals over the full range of the predicted values, which would suggest compliance with the assumption of homogeneity of variance (Figure S5).As an example, Figure 1 shows the sensitivity analysis for Scots pine and Maritime pine, respectively, in regard to the effect of the selected variables on tree height.Ho provided the largest contribution to the h-d models for both species, and showed that the tree heights increased with the increasing effect of stand development described by Ho and dq (Figure 1a,b,d,e).Moreover, the mixing proportion showed a different sign between species, indicating apparent differences of the effects of species mixing in the h-d fitted models (Table 2).When all of the other conditions were assumed to be the same, the tree height increased with increasing m PS for Scots pine (Figure 1c) and decreased with increasing m PT for Maritime pine (Figure 1f).Thus, Scots pine trees growing in monospecific stands reached larger heights than trees in mixed stands at a specific stand condition (including stand density).In contrast, an increase in m PT implies a smaller tree height for Maritime pine, so the heights associated with a given d were greater in mixed stands (Figure 2).effects of species mixing in the h-d fitted models (Table 2).When all of the other conditions were assumed to be the same, the tree height increased with increasing for Scots pine (Figure 1c) and decreased with increasing for Maritime pine (Figure 1f).Thus, Scots pine trees growing in monospecific stands reached larger heights than trees in mixed stands at a specific stand condition (including stand density).In contrast, an increase in implies a smaller tree height for Maritime pine, so the heights associated with a given d were greater in mixed stands (Figure 2).

Height-Diameter Function Validation
The most parsimonious h-d functions (Table 2) were validated based on the independent dataset, a set of pure and mixed plots selected from the third NFI plot network (Figure S3).Validation statistics were calculated and compared with the performance of the equations proposed by Lizarralde et al. [25,26] for monospecific stands of both species.When the fitted models were applied to trees in mixed stands, the Scots pine h-d function had a larger model efficiency and lower RMSE than of those of Maritime pine (Table 3).

Height-Diameter Function Validation
The most parsimonious h-d functions (Table 2) were validated based on the independent dataset, a set of pure and mixed plots selected from the third NFI plot network (Figure S3).Validation statistics were calculated and compared with the performance of the equations proposed by Lizarralde et al. [25,26] for monospecific stands of both species.When the fitted models were applied to trees in mixed stands, the Scots pine h-d function had a larger model efficiency and lower RMSE than of those of Maritime pine (Table 3).Results were consistent across stand composition when prediction bias was compared between the models fitted in this study and those proposed by Lizarralde et al. [25,26].For both species, the fitted model that included mixed-species stands attributes reduced the total prediction error in the mixed stands, increasing the model efficiency by 6% and 1%, for Scots pine and Maritime pine, respectively.The gain in efficiency in pure stands was about 2% for both species.The fitted models reduced the mean bias of tree height predictions for both species in mixed and monospecific stands across medium to small diameter classes (d < 25 cm; Figure 3).Bias from both models in other diameter classes followed the same patterns.Bias was evident and sizable in larger diameter classes (d > 65 cm) for Maritime pine, although these diameter classes were outside the range of the values that were used to fit the models (Table 1).diameter classes followed the same patterns.Bias was evident and sizable in larger diameter classes (d > 65 cm) for Maritime pine, although these diameter classes were outside the range of the values that were used to fit the models (Table 1).

Basal Area Increment Model
Differences in growth rates between mixed and pure stands were more evident in Maritime pine that in Scots pine (Figure S6).In the fitting model, the inclusion of random effects greatly improved the goodness-of-fit indicators, decreasing AIC and producing significant differences in the likelihood ratio.For both species, a power variance function depending on d performed well for homogenizing the variance.A visual check of the model residuals indicated no major departures from the assumptions of normality and homogeneity of variance.
According to the base model structure, the basal area growth increased with the initial tree diameter, d, and with CR, but the magnitudes differed between species.The site index did not improve the performance of the models in either species, as indicated by non-significant differences from the likelihood ratio tests, so it was not included in the base model.Rh had high collinearity with the diameter (variance inflation factor, VIF > 10), so it was also excluded.The inclusion of competition in the model structure considerably improved the BAI models for both species.

Basal Area Increment Model
Differences in growth rates between mixed and pure stands were more evident in Maritime pine that in Scots pine (Figure S6).In the fitting model, the inclusion of random effects greatly improved the goodness-of-fit indicators, decreasing AIC and producing significant differences in the likelihood ratio.For both species, a power variance function depending on d performed well for homogenizing the variance.A visual check of the model residuals indicated no major departures from the assumptions of normality and homogeneity of variance.
According to the base model structure, the basal area growth increased with the initial tree diameter, d, and with CR, but the magnitudes differed between species.The site index did not improve the performance of the models in either species, as indicated by non-significant differences from the likelihood ratio tests, so it was not included in the base model.Rh had high collinearity with the diameter (variance inflation factor, VIF > 10), so it was also excluded.The inclusion of competition in the model structure considerably improved the BAI models for both species.Competition terms relying on BA and BAL yielded the most parsimonious models compared to competition expressed in terms of relative stand density index.
Table 4 compares the different models that included symmetric or asymmetric competition structures.Not all the possible combinations of competition structure are listed, because non-significant parameters were excluded from the models, and the remaining parameters were re-estimated.Clear differences between species in terms of competition structure were observed after ranking models by AICc and Akaike weights (w i ).For Scots pine, the best ranked model included both size-symmetric and size-asymmetric for predicting the basal area tree growth in mixed and pure stands, with species composition adding little predictive ability (Table 4).The model containing species-specific size-asymmetric competition, BAL intra and BAL inter , increased AICc very little (∆AICc = 1.77), endorsing the idea that species mixing may not have a strong influence on Scots pine basal area growth.The results for Scots pine contrasted with the best basal area periodic annual increment (cm 2 5 years −1 ) (BAI) model for Maritime pine.Here, the effects of competition on the basal area growth rate in the most parsimonious model were most accurately predicted by the combination of size-symmetric and size-asymmetric indices separated into intraspecific and interspecific components (Table 4).Akaike weights, w i , for the model separating asymmetric competition into species-specific components indicated the relative superiority of the second-ranked model containing BAL with no discrimination between competitor species, and was much better than the model with no separation of total competition between species (BA or BAL).

Table 4. Ranking of basal area increment models comparing alternative measures of competition.
Competitive status was compared for the size-symmetric (BA) and size-asymmetric competitors (BAL) separated into intraspecific and interspecific components.Models ranked according to their ∆AICc: difference of AICc values between the best model and the ith model; and w i : Akaike weights.The BAI models with the lowest AICc value were fitted by the REML procedure (Table 5).All of the parameters were significant at p-value < 0.05.Species competition relationships only affected BAI in Maritime pine.Here, the symmetric competition effect was imposed mainly by other Maritime pine, and the symmetric interspecific competition imposed by Scots pine was not significant.Thus, at increasing proportions of basal area of Scots pine (decreasing BA intra ) in mixed stands, Maritime pine BAI increased (Figure 4a).On the other hand, the size-asymmetric competition of Scots pine was greater than the intraspecific competition, implying an increased competition effect when the larger competitors were Scots pine rather than Maritime pine competing against itself.This relationship resulted in a decrease of Maritime pine BAI with an increasing BAL proportion of Scots pine in a mixed stand (Figure 4b).Species mixing imposed a positive effect in size-symmetric competition, but a negative effect on size-asymmetric competition.

Species
Table 5.Estimated parameters (standard errors in brackets) and fitting statistics for the final basal area increment models of the two species.Parameters were estimated by the restricted maximum likelihood method (REML).Variables are defined as follows:σ 2 ε , residual variance; σ 2 i , variance for the plot random effects; δ, parameters of the variance function (Equation ( 5)).Bias, EF, and RMSE are the mean bias model efficiency and root mean square (Equations ( 8)-( 10)), respectively.

Scots Pine Maritime Pine
Fixed parameters the plot random effects; , parameters of the variance function (Equation ( 5)).Bias, EF, and RMSE are the mean bias model efficiency and root mean square (equations ( 8)-( 10)), respectively.

Scots Pine Maritime Pine
Fixed parameters   Figure 5 depicts the effects of species composition of both size-symmetric and size-asymmetric competition on the BAI of Maritime pine, comparing different scenarios when trees were influenced by intraspecific competition only, pure stands (Proportion = 1), and when competition was the combination of intraspecific and interspecific competition, mixed stands (Proportion < 1).The positive effect species mixing on Maritime pine tree growth was smaller as the proportion of Scots pine in the size-symmetric competition increased.However, BAI reduced strongly as the total size-asymmetric competition of Scots pine increased (e.g., in codominant and suppressed trees).The Scots pine showed no species-mixing effects, and the influence of BA, BAL, and CR on the BAI model was simulated in Figure S7.

Discussion
Contrasting mixing effects on two or more species in mixed forest stands are common among a range of forest types [9,10,13,50,51].Our results similarly showed that species-mixing had contrasting effects between species in both the generalized h-d and basal area growth models.Both models described the stand structure and stand dynamics of mixed stands relying on temporary plots grouped in triplets and located where two species such as Scots pine and Maritime pine mix, as is the case in the Northern Iberian Range.Mixed and pure plots covered the prevalent stand development stages, site conditions, and mixing proportions in the region, and were exposed to the same silvicultural regime within a given triplet.Typically, a triplet design allows the control of some stand conditions analogous to blocking, avoiding confounding factors at least for part of the stand development [52].In this way, estimates of tree height and BAI would accurately represent the species interactions in models that are sufficiently flexible to incorporate mixing effects and test for their magnitude and direction.

Influence on Species-Mixing in Tree Allometry
For both species, h-d models described a large proportion of the variation in tree height among plots and exhibited no substantial trends in the residuals to indicate bias (Figures S4 and S5).This result suggested that the selected model (Equation ( 4)) was sufficiently flexible to accommodate mixed species stands and mixing effects that influence the basic structural relationships.Although the final fitted models for mixed species stands differed from the basic formulation of the models developed by Lizarralde et al. [25,26] across a broader geographic area, both included stand variables that reduced bias and increased the precision of the model estimates [38,41].Testing predictions from these new equations on an independent dataset showed modest increases in the accuracy of height estimates for trees growing in mixed and pure stands (Table 3 and Figure 2).Dominant height was important to describe the effect of site quality as defined by stand growth and yield capacity [38,39], and it was the most influential covariate in the h-d models (Figure 1a,d).All else being equal, dq measured the degree of inter-tree competition, and secondarily the stage of stand development [41,53,54].Although the same variables have already been used in several h-d functions for multi-species forest [39,40,50,54], in those studies, the effects of species interactions on h-d relationships were not quantified.Here, we show that for mixed stands, the explicit evaluation of mixing proportion based on species-specific growing space requirements could be advantageous for modeling the more complex effects of the vertical structure of mixed species stands compared to pure stands.Different species almost certainly differ in their influence on competition for light in a forest canopy given their variation in potential crown dimensions, branching structure, leaf area density, and light capture [17,55,56].Although other tree and stand variables such as stand age or spatially explicit competition indices may also help fine-tune h-d relationships [39], these were not considered because they require a substantial increase in the cost of forest inventory, which restricts the practical application of the developed predictive functions.
The structure of Equation ( 4) considers Ho and m as complementary parameters, with Ho acting as an asymptote, and m quantifying relative species prevalence in the h-d models.The effect of m varied between the species in the fitted models (Table 2), and the sign of their parameters and the size ranges of the two species in mixed species stands (Table 1) seemed to support the tendency of Maritime pine to overtop Scots pine in mixed stands (Figure 2 and Figure S8).The varying magnitude of the effects of species-mixing in the h-d relationship observed for other mixtures [55,57,58] suggest that the degree of differences in the h-d relationship might depend on shade tolerance, crown width potential, relative juvenile and mature height growth patterns, and maximum potential stand density, with a greater response exhibited by less shade-tolerant species [9].
The predictions from the h-d models fitted to the triplets data for both species perform better than estimates from the previously available h-d models for pure stands [25,26].Besides the inclusion of terms that expressed the species-mixing effect on the h-d models, consideration of the hierarchical structure of the data in the NLME modeling approach increased the proportion of variation in height, which was accounted for by including random plot effects in addition to only fixed-effects nonlinear models [38,42].Increasing the accuracy of tree height estimates in mixed stands has important implications, because differences in tree morphology directly affect crown competition, stem volume, biomass production, mechanical stability, and wood quality predictions.In addition, small differences in initial stand structure associated with differences in species composition and proportions of the species are magnified over time, and impose profound changes in stand dynamics [59].
In previous analyses, mixed stands of Scots pine and Maritime pine showed changes in their vertical structure compared with pure stands, which resulted from species interaction rather than the simple combination of species with different traits and structural morphology [33].The heterogeneity in the vertical structural of mixed stands has been related to over-yielding at the stand level, and sometimes greater stand productivity than in single-species stands [8,14].These results underscore the need to develop models that reliably account for species interactions as part of the stand dynamics of mixed-species stands.

Tree Basal Area Growth Models Adapted to Mixed Stands
Our analysis showed that interspecific competition was only relevant for Maritime pine BAI; however, these effects differed depending on whether the competition was symmetric or asymmetric.In contrast, Scots pine BAI was unaffected by mixing, i.e., insensitive to species composition of BA.
For Maritime pine, only intraspecific size-symmetric competition was significant, which might suggest that an increase in BAI can be achieved by a reduction in competition for the below-ground resources in mixtures with Scots pine compared to homogeneous pure stands.However, this positive effect might be overridden by the greater negative effect of interspecific size-asymmetric competition, which might indicate stronger competition for light [46].Figure 5 shows the effects of interspecific competition and competitive reduction occurring simultaneously when comparing size-asymmetric and size-symmetric competition depending on the competitor species.This simulation also indicates that species exhibited unbalanced competitive relationships for above-ground and below-ground resources, resulting in a positive species-mixing effect for Maritime pine.These results are consistent with the over-yielding observed at the stand level analyzing the same triplets [33], and with the species competitive relationships observed for Maritime pine using NFI data [34].The different spatial scale of the data used in both analyses might explain these discrepancies.
Although deciphering the underlying ecological causes of these interactions were not the aim of this study, we can propose some possible explanations.The species-specific plastic responses of physiological and morphological traits to shade conditions could be one driving mechanism [60,61].Differences in shade tolerance between both species [62] might facilitate the complementary crown architectures [63] and canopy vertical stratification [33] observed in mixed stands of Maritime pine and Scots pine.Thus, intercepted light can also be used more efficiently, suggesting that light-related interactions may contribute to the species-mixing effect on stand productivity.However, in mixed stands of light-demanding species such as pines, the degree of differential preference for light is probably not the single controlling mechanism of the positive [64], negative [65], or null [66] effects observed in growth.Complementarity in traits for securing below-ground resources such as differences in water-use efficiency [67] or varying growth responses to fluctuating or directional drought conditions [68] needs special attention.
That the age dependent site index was not significant in any final basal area growth models probably indicates the absorption of site quality into random plot effects.The unimportance of site index in this analysis also coincided with the results of stand-level models developed for the same species in Spain [28,69].Initial basal area may be enough to indicate stand quality, and might produce accurate productivity estimates, especially in stands exposed to low levels of silvicultural intervention [28], such as pure and mixed stands located in the Northern Iberian Range.Furthermore, species-specific site index equations based on stand age are usually not available in multi-species stands [16], because interspecific competition probably suppresses top height growth for different stages in stand development [17,33], altering the height-age relationship and biasing SI estimations [56,70].Although a site index based on reference diameter is available in Spain [71], caution is advised in its application because of the inherent sensitivity of diameter growth to stand density relative to height growth [72,73], especially in mixed-species stands where stand density can alter growth dynamics [66,74].A topic for further exploration might also be compensatory growth in one species or the other along gradients in productivity that parallel changes in species composition and fluctuations in annual growing conditions.

Implication for Tree-Level Model in Mixed Stands
The generalized h-d and individual-tree BAI models developed in this study can be used for Scots pine and Maritime pine trees growing in pure and mixed stands.The effect of the interspecific competitive environment was successfully represented in the structure of the models.Combining the NLME modeling approach with parameters expressing species mixing enhanced the performance of tree height and basal area growth estimates compared with the previous models fitted for pure stands.The superior performance of multilevel mixed-effects models over fixed-effects models fitted to data collected under a more complex sampling design has already been widely discussed in the literature [38,42].However, the magnitude of variations in allometric relationships and tree growth due to species-mixing is a subject of intense debate, with evidence of adverse and positive effects that depend on multiple factors operating simultaneously [17,64,75,76].Therefore, the application of these models is restricted to the range of data as well as to the Northern Iberian Range, because species-mixing effects can vary among regions [33,34].Climate-sensitive models could have broader applications for h-d relationships and tree growth functions [77,78]; however, analyzing climatic effects is beyond the scope of this study, both because of the narrow geographical amplitude of the triplet network and the relatively short growth series.However, further work is needed to test if the specific years and annual fluctuations in climatic variables might affect the basal area growth and the h-d relationships, or shift the effect of species-mixing on productivity and tree allometry.
The adaptation of IBERO model to mixed stands can be used as a decision-support tool for evaluating the impact of different silvicultural options in the context of sustainable forest management.In this study, we proposed new species-specific models for adapting h-d functions and growth sub-modules of the tree-level IBERO model to mixed stands of Scots pine and Maritime pine in the Northern Iberian Range.Both h-d relationships and basal area growth equations might have larger effects on tree-level models in terms of gains or losses in productivity due to species interactions [33,34,63].Adapting growth models developed for monospecific stands to mixed stands requires that they become sensitive to species composition [2], and capture sufficient accuracy in upscaling species-mixing effects from the tree-level to the stand-level.For instance, volume and volume growth estimates in mixed stands that rely on volume equations developed for monospecific stands might result in overestimation or underestimation of productivity if form differs substantially [9,16].
Obviously, other IBERO sub-modules require analysis to determine the species-mixing effects on the complete model structure and evaluate the performance of growth estimates on mixed stand simulations.Several studies have demonstrated that species-mixing can have significant effects on species crown size and shape [63,79], taper equations [9], height growth [56], and mortality rates [11], all of which could influence the productivity of mixed-species stands.A wider range of stand structures is clearly possible in mixed-species stands, creating a challenge to both sampling strategies and the choice of model form for individual-tree equations.Therefore, conversely, caution must be exercised that mixing effects are not artefacts of inappropriate model structure or the biased quantification of mixing proportions.The IBERO site quality sub-module interacts with different components within the model architecture and relies on age-related site index curves developed for both species in monospecific conditions.Although the side index was not included in the formulation of the model that was developed in this study, several authors have demonstrated that species interaction effects can change in response to site conditions [15,78].As mentioned before, using the site index for even-aged monospecific stands could be ineffective in mixed-species stands.Hence, alternative methods to evaluate site quality based on geocentric measures [70,80,81] could be used to establish an alternative productivity index for mixed-species stands, especially in stratified mixtures in which site index curves have not been developed for the more dominant species, which is the least likely to have been overtopped throughout stand development.

Conclusions
Compared with previous studies, the originality of our work lies in that our h-d and basal area growth models explicitly accounted for species mixing effects in addition to tree and stand variables that describe the tree competitive environment and stage of stand development.The integration of variables reflecting the species-mixing effects in models to estimate tree height and basal area increment in mixed stands were shown to enhance the performance of predictions compared to available models parameterized for Scots pine and Maritime pine in pure stands.Both models could be incorporated into the tree-level model IBERO for the estimation of stand dynamics, growth, and yield of mixed stands.

Figure 1 .Figure 1 .
Figure 1.Effects of dominant height (Ho), quadratic mean diameter (dq), and species mixing proportion of target species ( or ) on the final h-d models in Equation (4) for Scots pine (a-c) Figure 1.Effects of dominant height (Ho), quadratic mean diameter (dq), and species mixing proportion of target species (m PS or m PT ) on the final h-d models in Equation (4) for Scots pine (a-c) and Maritime pine (d-f).Curves were produced using the parameter estimates in Table 2 and varying one explanatory variable at a time.Mean values of the data were used for predictors and the range of the variable of interest in the figure.

Forests 2018, 9 ,
x FOR PEER REVIEW 10 of 21 and Maritime pine (d-f).Curves were produced using the parameter estimates in Table2and varying one explanatory variable at a time.Mean values of the data were used for predictors and the range of the variable of interest in the figure.

Figure 3 .
Figure 3. Mean bias (m) by 5-cm diameter class from the best generalized h-d models fitted to the triplet data (red circles and lines) and from the Lizarralde et al. [25,26] h-d functions (blue circles and lines).

Figure 3 .
Figure 3. Mean bias (m) by 5-cm diameter class from the best generalized h-d models fitted to the triplet data (red circles and lines) and from the Lizarralde et al. [25,26] h-d functions (blue circles and lines).

Table 1 .
Tree and plot-level attributes in mixed and pure plots for each species across triplets.

Table 2 .
Parameter estimates, variance components, and goodness-of-fit statistics for the final generalized h-d model selected for each species.All of the parameters were significant with p < 0.05, standard errors in brackets.∆AIC is the reduction in AIC achieved by including species mixing proportion in Equation (4).

Table 3 .
[25,26]ion statistics, based on data from National Forest Inventory (NFI) plots in the North Iberian Range, comparing the best generalized h-d models fitted to the triplet data and the equations proposed by Lizarralde et al.[25,26].

Table 3 .
[25,26]ion statistics, based on data from National Forest Inventory (NFI) plots in the North Iberian Range, comparing the best generalized h-d models fitted to the triplet data and the equations proposed by Lizarralde et al.[25,26].: root mean square error (m), EF: model efficiency; n: number of trees selected for the validation, * Lizarralde et al.[25,26]. RMSE