Productivity Estimations for Monospecific and Mixed Pine Forests along the Iberian Peninsula Aridity Gradient

National Forest Inventories (NFIs) are the primary source of information to fulfill international requirements, such as growing stock volume. However, NFI cycles may be “out of phase” in terms of the information required, so prediction techniques are needed. To disentangle the effects of climate and competition on stand productivity and to estimate the volume of stocks at national scale, it is important to recognize that growth and competition are species-specific and vary along climatic gradients. In this study, we estimate the productivity of five pine species (Pinus sylvestris, Pinus pinea, Pinus halepensis, Pinus nigra and Pinus pinaster), growing in monospecific stands or in mixtures along an aridity gradient in the Iberian Peninsula, based on Spanish NFI data. We study the stand volume growth efficiency (VGE), since it allows the comparison of volume growth in monospecific and mixed stands. The results reveal the importance of considering the aridity when assessing VGE. Moreover, it was found that, in general, admixture among pine species leads to modifications in the VGE, which can vary from negative to positive effects depending on species composition, and that this is always influenced by the aridity. Finally, we provide simple growth efficiency models for the studied pines species which are valid for both monospecific and mixed stands along the aridity gradient of the Iberian Peninsula.


Introduction
Spain, as a European country, must fulfill international requirements, which demand information about forest indicators such as those requested by the FAO (Food and Agriculture Organization of the United Nations) or UNFCC (United Nations Framework Convention on Climate Change), among others. One such requirement is the stock volume [1,2], which is required at five-to ten-year intervals, or in some cases annually [3]. These requirements are used in global reports, which serve as a basis for international policy making. Consequently, the NFIs are the primary source of information used for producing reliable forest-relate estimates [1]; however, due the different designs of the NFIs they must be adapted so that these international requirements can be met [4]. In Spain, as occurs in other NFIs, the cycles are longer than desirable for the purposes of international reporting, with each plot being re-measured every ten years. Moreover, although in other countries the plots measured each year cover all the national forest area, which allows the use of model-based estimation or similar techniques [5,6], in Spain the plots do not cover all the forest area, but rather, are located each year in a different specific region, making it necessary to update stock volume estimates in order to provide annual information at national level [7]. Hence, one simple option could be to estimate annual volume growth rates based on the growth observed in previous SNFI cycles. However, this method may be problematic due to the continuous changes in forest area over time, meaning that past volume growth information is not available for all plots measured in the current SNFI cycle. An alternative method is the use of growth models based on past data from SNFI, which might enable the estimation of annual volume growth for all SNFI plots, including new plots established in the last inventory.
Although several forest growth models have been developed using Spanish NFI data, most of them focus on tree growth, often for specific species or regions. However, for national-scale data, hardly any models have been developed at stand level [8,9]. Gómez-Aparicio et al. [10] modelled tree growth of the main Iberian forest species, revealing high variation in tree growth and sensitivity to competition along climatic gradients. These authors also pointed to the importance of identifying the specific species present with regard to competition effects. However, such variation in the patterns found at tree level might not always be relevant at stand level due to emergent properties [11]. Therefore, to estimate the volume of stocks at national level it is important to disentangle the variation in the patterns of climate and competition effects on productivity at stand level.
Forest productivity depends on environmental conditions, stand age or stage of development as well as on stand structure, particularly stand density and species composition [12][13][14]. To assess the effect of environmental conditions on productivity, the concept of site index has traditionally been used in even-aged monospecific stands [15]. However, this index is difficult to apply for more complex structures [16,17] or in those forest inventories, such as SNFI, where information on age is lacking. In these cases, growth models often include climate and soil variables to reflect the productivity variation with environmental conditions [16].
Stand density is a critical factor for forest productivity because stand growth usually increases with density before levelling off at higher densities [13,18,19]. Moreover, density can modify the effects of species mixing on productivity [20][21][22]. Hence, density is key in the assessment of both stand growth and species competition. Regarding the influence of species composition, numerous studies report a positive effect of species mixing on productivity or overyielding [23,24], although this effect depends on species composition and might vary with environmental conditions [25][26][27][28].
In this study, we focus on the productivity of the main pine species (Pinus sylvestris L., Pinus pinea L., Pinus halepensis Mill., Pinus nigra Arn. and Pinus pinaster Ait.) growing in monospecific and mixed stands across the Iberian Peninsula given their economic, social and ecological importance [29]. According to Montero and Serrada [30], pine species are distributed over more than five million ha as monospecific stands, and almost half a million ha in pine-pine mixtures. In terms of volume, the five target species analyzed produce about 47% of the total timber volume in Spain.
Focusing on the Iberian Peninsula, one of the most evident environmental gradients is that of aridity [31], which has been found to influence forest growth and mortality [10,32,33], as well as the maximum size-density relationship (MSDR), of the main pine species [34]. As regards the mixing effects on pine-pine mixtures, it might be expected that species mixing would not result in overyielding due to the similar traits of the species, which might limit the presence of niche complementarity. However, some studies point to overyielding in mixed Mediterranean pinewoods [35][36][37].
In this paper, the main objective is to provide equations to update pine VGE when data are not available (i.e., between SNFI cycles). For this, we aim to identify how the spatial variation in aridity along the Iberian Peninsula modifies the species productivity and whether it influences on the mixing effects. We hypothesized that (i) productivity decreases with site aridity, with variation in stand productivity along an aridity gradient being similar to the MSDR pattern, as higher carrying capacity is linked to greater productivity; (ii) this variation in productivity differs among species; and (iii) species mixing modifies stand productivity, with mixing effects being influenced by aridity. To estimate productivity, we used the stand volume growth efficiency (VGE) as defined by [22,38], which allows volume growth in monospecific and mixed stands to be compared.
The specific objectives were, therefore, to develop specific pine VGE models using data from the SNFI and to assess the effect of the aridity as well as the mixture of pine species on stand volume growth efficiency.

Data
The data used for this study came from two consecutives cycles of the SNFI for the whole Iberian Peninsula; the Second and Third NFI (SNFI-2 and SNFI-3) for all provinces except for Navarra, Asturias and Cantabria where, because of the lack of continuity of the information, the SNFI-3 and Fourth NFI (SNFI-4) were used. These inventories were carried out using permanent sample plots located at the nodes of a 1 km square grid, which were re-measured approximately every 10 years (SNFI-2 from 1986 to 1996, SNFI-3 from 1997 to 2007 and SNFI-4 from 2007 to the present). The SNFI permanent plots are composed of four concentric subplots with radii of 5, 10, 15 and 25 m, where diameters and heights of all the trees with breast-height diameters of at least 7.5, 12.5, 22.5 and 42.5 cm, respectively, were measured. Using the appropriate expansion factors for the concentric subplots, main stand variables (total and per species) were estimated from individual trees. Stand volume growth was estimated using the Martin [39] method. This method only requires the volume of the trees measured at any inventory to be known, regardless of the plot measurement radius. Although this produces some negative stand volume growth values due the changes in the expansion factor of growing trees, it ensures additivity, i.e., the difference between time-separated volume estimates is equal to the sum of growth component estimates.
The species studied were five pines native to the Iberian Peninsula (Pinus sylvestris L. (Ps), Pinus pinea L. (Pp), Pinus halepensis Mill. (Ph), Pinus nigra Arn. (Pn) and Pinus pinaster Ait. (Pt)). Plots located in monospecific stands of pine and mixtures of pine species were used in this study. Plots affected by silvicultural fellings in more than 5% of the total basal area were not considered, nor were plots growing at very low densities (Relative Density (RD) lower than 0.2, see Section 2.3. for further explanation of the RD concept). The presence of species different from the five target species was also limited (RD of non-pine species < 0.05). Species proportion by area (defined in Section 2.3.) was used as the criterion for species composition classification, i.e., monospecific vs mixed stand. Species origin (plantation or natural) was also considered using raster maps for each target species [40]. These maps are based on the Spanish Forest Map at 1:50,000 scale (https://sites.google.com/site/sigforestspecies/home/mapas-de-especies). Table 1 shows a summary statistic of the monospecific and mixed plots used in the study. N is the number of trees per hectare; dg is the quadratic mean diameter (cm); Ho is the dominant height (m); IV is the volume increment (m 3 ·ha −1 ·year −1 ), M is the Martonne aridity index (mm· • C −1 ); RD is the relative density. The subscript sp refers to the target species (in mixed stands those named first), while j refers to companion species (the second species named in the mixture). When a variable does not have any subscript, it refers to the plot. The number of plots is given under the species name (#).
Furthermore, mean annual temperature (T, • C) and annual precipitation (P, mm) were used to assess the aridity gradient. Both variables were obtained for each SNFI plot from raster maps with a 1 km resolution developed by Gonzalo Jimenez [41]. As a measure of the aridity, the Martonne aridity index [42], calculated as M = P/(T + 10) ( • C/mm), was used.

Volume Growth Efficiency Models
With the aim of studying the variation in productivity along the aridity gradient of the Iberian Peninsula, volume growth efficiency (VGE) models per pine species (sp) were developed. The volume growth efficiency (VGE) was estimated by the volume growth of the species sp, IV sp (m 3 ·ha −1 ·year −1 ), divided by the proportion of this species sp (p sp ), (see Section 2.3. for species proportion definition). This variable was preferred as a response variable in order to accommodate the mixing effects among pine species, as it permits the comparison of species growth in monospecific and mixed stands [38].
In a first step, we developed general basic VGE models per species, linearized using natural logarithm but without considering species composition (Equation (1)). To avoid bias in the models, the negative values of IV sp , which appears due the concentric sample plot design, were also included by correcting the VGE values as VGE' sp = VGE sp + 1 − min (VGE sp ). ln(VGE' sp ) = a 0 + a 1 · ln(dg sp ) + a 2 · ln(Ho) + a 3 · ln(RD) + a 4 · ln(M) + a 5 · Origin (1) where for each plot, dg sp was the quadratic mean diameter of the species sp (cm); Ho, the dominant height (m); RD the relative density; M the Martonne aridity index (mm • C −1 ); and Origin was a dummy variable, the value of which was 1 when the species sp was a plantation, and 0 when it was a natural species. In these models, dg sp and Ho surrogated the stand age or development stage and the site index, respectively [22,43]. The min (VGE sp ) for all species are in Table S3.
In a second step, in an attempt to determine whether there are significant mixing effects that cause under or overyielding in productivity, the basic model was modified considering the presence of the other pine species (j) different from the target species. The variables used to assess the mixing effects were the relative densities of the pine companion species (RD j ). Additionally, interactions of these variables with M were included in order to test whether the mixing effects on VGE were influenced by aridity, in accordance with the approach by [28]. This total model, considering mixture effects, is presented in Equation (2). ln(VGE' sp ) = a 0 + a 1 · ln(dg sp ) + a 2 · ln(Ho) + a 3 · ln(RD) + a 4 · ln(M) + a 5 · Origin + a 6j · RD j + a 7j · M · RD j (2) Equation (2) represents a single model for each target species sp, and is suitable for both monospecific stands and for every possible mixture of the target species with the rest of the pine species. Thus, for each companion species j, RD j represents the relative density of this species in the plot, its value being equal to zero in monospecific stands. For each target species, the correlation between independent variables included in the model has been analyzed.

Relative Density and Species Proportions by Area
The species Relative Density (RD sp ) was calculated for each plot as the ratio between the current (N sp ) and the maximum (N max_sp ) number of trees per hectare (Equation (3)).
log(N max_Pn ) = 11.819 + 0.209 · log(M) − 1.8759 · log(dg) (7) log(N max_Pt ) = 9.784 + 0.937 · log(M) The total stand relative density (RD) per plot was calculated as the sum of the RD sp of all the pine species present in the plot, i.e., RD = RD sp . The RD and RD sp were used to estimate species proportion by area (p sp ) [17]: This definition allows mixtures to be analyzed, otherwise miscalculations could occur in the comparison between monospecific and mixed stands [44,45].
For each species sp, the species proportion by area p sp was used to calculate VGE sp = IV sp /p sp and to classify plots as monospecific, when the proportion of other pine species is null, p j = 0, or as pine-pine mixtures when the proportions of the two pine species are different from zero.

Model Fitting
Due to the fact that during the SNFI there were several teams doing the field work in different Spanish provinces, as well as in different years or seasons, mixed models were used to consider possible between-plot correlation. A random term was included in the intercept with the province as grouping structure.
The basic models (Equation (1)) and models including mixing effects (Equation (2)) were fitted using the lme procedure in the nlme package [46], for the R software [47]. The coefficients were included in the model only if they were statistically significant (p value < 0.05). Variables describing mixing effects and their interactions with M, were incorporated into the models if they led to a statistically significant improvement in the quality of the model fit as assessed at the alfa = 0.05 significance level using the F-test based on the extra sum of squares principle.
Comparisons between models according to Equation (1) and Equation (2) were performed by calculating the Akaike Information Criteria (AIC), while marginal and conditional R 2 were used to determine goodness-of-fit [48]. The marginal R 2 describes the proportion of variance explained by the fixed factor(s) alone, while conditional R 2 describes the proportion of variance explained by both the fixed and random factors [49,50].
The bias induced by logarithmic transformation was corrected for predictions by using the Beauchamp and Olson [51] coefficient (Table S3).

Uncertainty Assessment
SNFI data do not provide the high level of control of the conditions (e.g., site or stand structure) necessary to ensure that the effects of mixture on productivity are detected. The analysis using the models developed from this dataset only provide statistical significance of variables, which makes it difficult to determine whether they are a consequence of mixture or a consequence of other hidden drivers, such as differences in the site, or different measurement periods for monocultures and mixtures. To address these deficiencies and, at least partially, to assess the consequence of not controlling all site conditions, a Monte Carlo bootstrap approach, following the methodology described by Condés et al. [28], was developed to evaluate the uncertainty associated with mixing effects.
For each pine-pine mixture (Table 1), a total number of 1000 replications were performed. In each replication a set of pseudo-triplets was selected (randomly and with replacement) which consisted of the same number of plots located in monospecific and in mixed stands. The model in Equation (2) was then refitted to the selected data. Predictions from this model were used to estimate mean values of VGE. The overall means of the 1000 replications and their corresponding standard errors were used to estimate the confidence interval of estimates at the 0.95 confidence level.

Results
As adding variables related to the mixing effect to the VGE models led to better goodness-of-fit, we only analyzed these models (i.e., Equation (2)) for the rest of the study. The models for all species improved using Equation (2), although the species which showed the least improvement were Ps and Pn. The results of Equation (1) (basic models without mixing effect), as well as the goodness-of-fit for the Equation (1) and Equation (2) model comparisons, can be found in Table S1 and Table S2, respectively. Table 2 shows the parameter estimates together with their standard errors for the final models according to Equation (2). Differences between marginal and conditional R 2 show the improvement resulting from including the province random effect, the greatest improvement found being that for Pn while the least improvement was found for Ph. The marginal R 2 values vary from 0.2196 to 0.4027, while the conditional R 2 varies from 0.2962 to 0.5220, the lowest values being obtained for Pp and the greatest values for Pt. Table 2. Estimated coefficients and standard errors for volume growth efficiency models (Equation (2)) for the five target species.  Figure 1 represents the variation in VGE for monospecific stands while all the independent variables remain constant except the one represented on the OX axis. In general, the highest VGE were found for Pt while the species with the lowest growing efficiency was Ph. The other three species, Pp, Ps and Pn, presented intermediate and similar efficiencies.
The quadratic mean diameter, as a surrogate of age or stage development, had a strong negative effect on VGE (Figure 1a), with the VGE of Ps being the most influenced by changes in dg, ( Table 2). Pt displayed a similar behavior while the species least affected by this variable was Ph, which is also characterized by having the shortest range of dg values. In general, the dominant height variation results in less steep VGE curves, especially for Pp, Ps and Pn. This effect was always positive, reflecting the site quality. Ph, despite the lower range of Ho values, showed a large variation in VGE (Figure 1b). The effect of density (RD) on VGE, although always positive, differed among species, with Ph being the species least affected by RD, although in general it was growing at lower densities ( Table 2).  The quadratic mean diameter, as a surrogate of age or stage development, had a strong negative effect on VGE (Figure 1a), with the VGE of Ps being the most influenced by changes in dg, ( Table 2). Pt displayed a similar behavior while the species least affected by this variable was Ph, which is also characterized by having the shortest range of dg values. In general, the dominant height variation results in less steep VGE curves, especially for Pp, Ps and Pn. This effect was always positive, reflecting the site quality. Ph, despite the lower range of Ho values, showed a large variation in VGE (Figure 1b). The effect of density (RD) on VGE, although always positive, differed among species, with Ph being the species least affected by RD, although in general it was growing at lower densities ( Table 2). Furthermore, the species origin also modifies VGE in all cases except for Pp (Table 2). When the species considered are plantations (Origin = 1), this variable multiplies the model by the estimated parameter (e ^ (a5 · 1)), which is the model multiplied by 1.18, 1.07, 1.17 and 1.11 for Ps, Ph, Pn and Pt, respectively.

Influence of Aridity on the Productivity
Of particular interest was the influence of the aridity on the growth efficiency of the different pine species growing in monospecific stands. It is important to note that, as hypothesized, the productivity of all pine species increased as the aridity decreased, i.e., higher values of VGE for higher values of the Martonne aridity index (M) (Figure 1d). The intensity of this influence was, however, much greater in the case of Pt than for the other four studied species, with Pt also being the species growing across the widest range of aridity conditions. In contrast, Ph is the species which is distributed along the shortest gradient of M values, although the influence of aridity on the productivity is similar to that of the other species (Table 2). Nevertheless, the aridity also has an indirect effect on productivity through the RD, as the maximum density also depends on the aridity. Figure 2 shows that the influence of RD on VGE is modified by aridity. The widest range of M values was for Ps, but the species that displays the greatest increase in VGE with an increase in M is Pt. Furthermore, the species origin also modifies VGE in all cases except for Pp (Table 2). When the species considered are plantations (Origin = 1), this variable multiplies the model by the estimated parameter (eˆ(a 5 · 1)), which is the model multiplied by 1.18, 1.07, 1.17 and 1.11 for Ps, Ph, Pn and Pt, respectively.

Influence of Aridity on the Productivity
Of particular interest was the influence of the aridity on the growth efficiency of the different pine species growing in monospecific stands. It is important to note that, as hypothesized, the productivity of all pine species increased as the aridity decreased, i.e., higher values of VGE for higher values of the Martonne aridity index (M) (Figure 1d). The intensity of this influence was, however, much greater in the case of Pt than for the other four studied species, with Pt also being the species growing across the widest range of aridity conditions. In contrast, Ph is the species which is distributed along the shortest gradient of M values, although the influence of aridity on the productivity is similar to that of the other species (Table 2). Nevertheless, the aridity also has an indirect effect on productivity through the RD, as the maximum density also depends on the aridity. Figure 2 shows that the influence of RD on VGE is modified by aridity. The widest range of M values was for Ps, but the species that displays the greatest increase in VGE with an increase in M is Pt.

Influence of Mixing Effects on the VGE
The effects of interspecific competition (RDj) on the VGE and the influence of the Martonne aridity index (M) differ among mixtures, with the effects not always being statistically significant ( Table 2). Figure 3 shows the effect of mixing on the VGE for each species, i.e., the comparison between VGE in mixed and monospecific stands of this species, and how it varies along the M aridity gradient along which the species is distributed. The magnitude of the effects differs among mixtures, being very small in some cases, such as Ps-Pn, even though the coefficients are significant (Table 2, Figure 3a). The mixing effect on productivity can be positive or negative depending on the mixture, but in all cases, the greater the proportion of companion species the more evident the effect.

Influence of Mixing Effects on the VGE
The effects of interspecific competition (RD j ) on the VGE and the influence of the Martonne aridity index (M) differ among mixtures, with the effects not always being statistically significant ( Table 2). Figure 3 shows the effect of mixing on the VGE for each species, i.e., the comparison between VGE in mixed and monospecific stands of this species, and how it varies along the M aridity gradient along which the species is distributed. The magnitude of the effects differs among mixtures, being very small in some cases, such as Ps-Pn, even though the coefficients are significant (Table 2, Figure 3a). The mixing effect on productivity can be positive or negative depending on the mixture, but in all cases, the greater the proportion of companion species the more evident the effect.

Influence of Mixing Effects on the VGE
The effects of interspecific competition (RDj) on the VGE and the influence of the Martonne aridity index (M) differ among mixtures, with the effects not always being statistically significant ( Table 2). Figure 3 shows the effect of mixing on the VGE for each species, i.e., the comparison between VGE in mixed and monospecific stands of this species, and how it varies along the M aridity gradient along which the species is distributed. The magnitude of the effects differs among mixtures, being very small in some cases, such as Ps-Pn, even though the coefficients are significant (Table 2, Figure 3a). The mixing effect on productivity can be positive or negative depending on the mixture, but in all cases, the greater the proportion of companion species the more evident the effect.  According to our results, mixing effects are always modulated by the aridity, i.e., for all species where RDj is significant, it appears interacting with M ( Table 2). In most of the mixtures, the humidity amplifies the mixing effect on productivity, making it more positive as for Pn-Pt, or more negative as for the admixtures Ps-Pn, Pp-Pt, and all admixtures where Ph is the companion species (Pp-Ph, Pn-Ph and Pt-Ph). In other cases, the mixing effects on productivity switch from positive to negative depending on the aridity, as for Ph-Pt and Pt-Ps. In these admixtures, the VGE of the target species is greater than in monospecific stands in the arid sites, while the opposite occurs when the humidity increases. According to our results, mixing effects are always modulated by the aridity, i.e., for all species where RD j is significant, it appears interacting with M ( Table 2). In most of the mixtures, the humidity amplifies the mixing effect on productivity, making it more positive as for Pn-Pt, or more negative as for the admixtures Ps-Pn, Pp-Pt, and all admixtures where Ph is the companion species (Pp-Ph, Pn-Ph and Pt-Ph). In other cases, the mixing effects on productivity switch from positive to negative depending on the aridity, as for Ph-Pt and Pt-Ps. In these admixtures, the VGE of the target species is greater than in monospecific stands in the arid sites, while the opposite occurs when the humidity increases.

Uncertainty Assessment
The above results point to the mean effects of mixture on productivity. However, although they are statistically significant, they could be the result of site differences between plots located in monospecific and mixed stands. The uncertainty analysis shows the standard errors for mixing effects on the VGE of each species. Figure 4 shows that the higher the RD the wider the standard errors of the mixing effect estimates, and the confidence intervals included the horizontal line equal to one (meaning no effects).
on the VGE of each species. Figure 4 shows that the higher the RD the wider the standard errors of the mixing effect estimates, and the confidence intervals included the horizontal line equal to one (meaning no effects).
The confidence intervals show that mixing effects differ from zero, i.e., confidence intervals do not include a horizontal line (Figure 4) for Pp when mixed with Ph, for Pn when mixed with Pt and for both species in the Pt-Ph mixtures. In general, although the confidence intervals are wider for the more humid conditions, i.e., higher M values, the above results are still valid. It is interesting, however, to highlight the influence of M on the effect of Pt on Ph species. When this admixture is located in more arid sites the overyielding in Ph is evident, but when humidity increases the effect of Pt on Ph is negative, i.e., underyielding, although with high uncertainty.
In the admixtures Pp-Pt, Pn-Ph and Pt-Ps the uncertainty analysis showed that although there are significant effects in the models, the confidence intervals indicate that in some cases, i.e., under particular conditions, the effect could be null.

Aridity-Dependent VGE Models
The volume growth efficiency models for the five pines species had a similar structure, differing only in the component expressing the species mixing effects (see Section 4.2). The variables used in VGE models have been successfully used in previous studies when analyzing stand productivity based on Spanish NFI data [22,28,37]. The effects of dominant height and quadratic mean diameter on VGE confirmed their suitability as surrogates of site productivity and stage of development, as reported in other studies [37,[52][53][54]. Considering both dg and Ho, the species least affected by both variables was Pp, probably due to the low stand densities promoted by forest management [55]. However, Ph was the species that presented the most extreme responses to both variables, that is, it was most affected by Ho, while it was least affected by dg. This reveals that despite its frugality [56], Ph benefits from an improvement in environmental conditions [57].
In our VGE models, the effects of site conditions were also considered through the inclusion of climatic variables [58,59]. We used the Martonne index for all species and admixtures as it is a simple index, which has given good results in previous studies [24,60], and a common variable which allows direct comparisons among species compositions. However, other studies used different climatic variables, such as evapotranspiration or mean temperature of the coldest or hottest month [8,61,62]. Hence, it is possible that the use of different climatic variables for each species would provide better results [10]. The magnitude of the effect of M is not completely linked to the M range of the distribution of the species. Pt has the largest distribution of M and is also the most affected by M. However, the distribution range of M in the case of Ps is similar to Pt, and the effect of M on Ps is almost the lowest. The greater variation in productivity along the aridity index for Pt is in accordance with the large geographical variation in tree growth and productivity [63][64][65][66][67].
Despite using two variables related with environmental conditions, M and Ho, the correlations between them, in all species, demonstrate that both indicate different aspects of local conditions. While M only considers climate conditions, Ho collects the site conditions as a whole (soil, elevation, aspect, etc.) ( Table S4). The effect of density on productivity has been extensively discussed in forest The confidence intervals show that mixing effects differ from zero, i.e., confidence intervals do not include a horizontal line (Figure 4) for Pp when mixed with Ph, for Pn when mixed with Pt and for both species in the Pt-Ph mixtures. In general, although the confidence intervals are wider for the more humid conditions, i.e., higher M values, the above results are still valid. It is interesting, however, to highlight the influence of M on the effect of Pt on Ph species. When this admixture is located in more arid sites the overyielding in Ph is evident, but when humidity increases the effect of Pt on Ph is negative, i.e., underyielding, although with high uncertainty.
In the admixtures Pp-Pt, Pn-Ph and Pt-Ps the uncertainty analysis showed that although there are significant effects in the models, the confidence intervals indicate that in some cases, i.e., under particular conditions, the effect could be null.

Aridity-Dependent VGE Models
The volume growth efficiency models for the five pines species had a similar structure, differing only in the component expressing the species mixing effects (see Section 4.2). The variables used in VGE models have been successfully used in previous studies when analyzing stand productivity based on Spanish NFI data [22,28,37]. The effects of dominant height and quadratic mean diameter on VGE confirmed their suitability as surrogates of site productivity and stage of development, as reported in other studies [37,[52][53][54]. Considering both dg and Ho, the species least affected by both variables was Pp, probably due to the low stand densities promoted by forest management [55]. However, Ph was the species that presented the most extreme responses to both variables, that is, it was most affected by Ho, while it was least affected by dg. This reveals that despite its frugality [56], Ph benefits from an improvement in environmental conditions [57].
In our VGE models, the effects of site conditions were also considered through the inclusion of climatic variables [58,59]. We used the Martonne index for all species and admixtures as it is a simple index, which has given good results in previous studies [24,60], and a common variable which allows direct comparisons among species compositions. However, other studies used different climatic variables, such as evapotranspiration or mean temperature of the coldest or hottest month [8,61,62]. Hence, it is possible that the use of different climatic variables for each species would provide better results [10]. The magnitude of the effect of M is not completely linked to the M range of the distribution of the species. Pt has the largest distribution of M and is also the most affected by M. However, the distribution range of M in the case of Ps is similar to Pt, and the effect of M on Ps is almost the lowest. The greater variation in productivity along the aridity index for Pt is in accordance with the large geographical variation in tree growth and productivity [63][64][65][66][67].
Despite using two variables related with environmental conditions, M and Ho, the correlations between them, in all species, demonstrate that both indicate different aspects of local conditions. While M only considers climate conditions, Ho collects the site conditions as a whole (soil, elevation, aspect, etc.) ( Table S4). The effect of density on productivity has been extensively discussed in forest science [12], with different response patterns depending on species, age, and site conditions [13]. In our models, all species presented a positive effect of relative density on VGE, the most affected species being Pn and Ps, and the least affected Ph ( Table 2). The magnitude of this effect is not related to their RD values, all species presenting similar values ( Table 1).
The RD was calculated using the maximum density lines depending on site conditions (M) developed by Aguirre et al. [34]. This means that the effect of RD on VGE was more evident in more humid sites (Figure 2), which agrees with results reported for P. sylvestris [68] and Picea abies (L.) [12]. The use of site-dependent species-specific maximum density is fundamental when the data used cover a large gradient of sites [34,69]. The density variation with site conditions modifies the competition equivalence coefficients and species proportions, which modifies the assessment of productivity [70].

Mixing Effects
The results show that, in general, when two pine species form a mixture, they have neutral or negative effects on each other (Figure 3), which can result in either similar productivity, or even underyielding in mixtures compared to monospecific stands (Table 3). Most previous studies of stand productivity in different forest mixtures point to positive mixing effects or overyielding [25,26,71], although under specific conditions underyielding also occurs [25]. Overyielding is linked to niche complementarity between species and/or facilitation, whereas underyielding indicates stronger interas opposed to intra-specific competition [72]. Our results indicate that competition is the main species interaction in pine mixtures, which may result from similarity in species traits [73]. However, for some mixtures and aridity conditions, positive effects of one species on the other were also found, as reported in other studies of Mediterranean mixed pinewoods [74,75]. Table 3. Summary of mixing effects on productivity at species level and resulting effects at stand level for low and high aridity (M). Net effects of mixing on stand productivity are influenced by environmental conditions, with contrasting effects depending on species composition [71]. With the exception of the Pt-Pn mixture, our results showed greater underyielding in more humid locations (Figure 3, Table 3). This pattern agrees with the stress-gradient hypothesis [76], which predicts that facilitation will occur in harsh conditions, and that it will be less pronounced in humid locations. In contrast, competition is the prevailing effect under more favorable conditions, that is, in more humid sites. A similar pattern with aridity was found for other mixtures of forest species in studies undertaken at tree level [33,77], as well in some studies conducted at stand level which used site index as an indicator of environmental conditions [25,44,54]. However, a recent meta-analysis reports greater overyielding with increasing humidity [26], in line with the results for the Pt-Pn mixture. These opposite patterns highlight the complexity of species interactions, which depend on the type of species and the availability of resources [78]. When a stand does not suffer from water deficit, the main limiting factor would be light, so overyielding can result from complementarity between species at aboveground level [20]. This may be the reason for the overyielding found in the Pn-Pt admixture. Few cases of complementarity or positive effects have been identified for the mixtures studied, which suggests low aboveground complementarity between species. This could be due to similar shade tolerances and crown structures [79], or to differences in growth patterns between species at a given site [80]. This may result in height dominance of one species over the other (explaining the neutral effect for one species (the dominant one) and the negative for the other).

Ps-Pn
Under less humid conditions the competition between the studied pine species is reduced, even displaying a certain degree of complementarity, i.e., positive effect ( Figure 3). For the Ph-Pt and Pt-Ps admixtures, the mixing effect changes from over-to underyielding under arid and humid conditions, respectively. This reflects the complexity of mixing effects, even for the same mixture [44], and the importance of considering site conditions when analyzing stand productivity [71].
Regarding the specific species compositions, our results differ partially from previous studies. Riofrío et al. [37] studied mixed stands of Ps and Pt and found that both species affect each other positively. In our study, however, we found that it was only Ps that affected the VGE of Pt, either positively or negatively depending on M. Trasobares et al. [81] also found that Ps and Pn affect each other at tree level, although the effect of Pn on the VGE of Ps was more evident, this being the only effect identified in our study for this mixture. When Pp and Ph interactions were studied by Cattaneo et al. [82], both species were found to affect each other negatively at tree level in a mixed plantation. Our results for the same admixture indicate that Ph affects the VGE of Pp, but not the other way around. This discrepancy in the findings may be due to differences between areas studied, since our analysis was conducted at national level, which implies a large variation in site conditions and this in turn will have an impact on the mixture effects [44]. Another difference that could explain the discrepancies is the scale of analysis, i.e., tree-or stand-level [83]. However, other methodological factors such as the type of data source [84] or the species proportion definition [44,85] could also influence the results.

Methodological Considerations
When using large-scale forest inventories to predict productivity, it is important to assess uncertainty [7,28,86], although it can be very sensitive to the data used [87]. In the present study, this technique has been used to analyze the effects of mixing from NFI data, as in Condés et al. [28]. NFI plots do not normally present high densities and consequently, the widest confidence intervals correspond to high densities for all admixtures (Figure 4), i.e., the least represented plots. Although for some mixtures the confidence intervals do not allow us to confirm significant species interactions, the analysis of uncertainty indicates the need to consider mixing effects when VGE is analyzed.
The suitability of NFI data for developing growth models has been questioned by some authors [27,35,[88][89][90]. The main critique is the lack of ceteris paribus control of conditions [27]. Therefore, the use of NFI data makes it difficult to determine whether the results obtained are due to mixing effects or to other factors which are beyond the scope of NFI data [28]. Another limitation is that NFI data does not allow differentiation between subspecies. In the cases of the species studied, Pn and Pt each have two subspecies (subspecies nigra and Salzmanii in the case of Pn, and atlantica and mesogeenesis in the case of Pt), the genetic variability of which could lead to different behavior [91,92]. In contrast, the main advantage of the NFI is that all species and admixtures are well represented with data systematically distributed along large gradients of both site and stand conditions [25,93,94]. This distribution of data is important to understand the effect of species mixing at larger scales [44].
Through the models obtained, the responses to international requirements in terms of volume can be more accurately provided for five of the main Spanish species. These models let to update growths, but to get more reliable updates it should take into account natural losses and fellings [95] using scenario analyses or past estimates of removals [95]. As data are collected in different years across the Iberian Peninsula, the use of these models allows us to extrapolate the results to a specific year at national level. However, the main limitation of the models developed is that they are only suitable for short periods of time, when both climatic conditions and management regimes do not vary [7,96]. Therefore, models are appropriated always when no unusual events, such as abnormal meteorological conditions, happens, which means additional uncertainties linked to climatic conditions under climate change [97].

Conclusions
The initial hypotheses were: (i) productivity decreases with aridity; (ii) this variation in productivity differs among species; and (iii) species mixing modifies stand productivity, with mixing effects being influenced by aridity.
Our models for estimating the productivity of the five main pine species (Pinus sylvestris L., Pinus pinea L., Pinus halepensis Mill., Pinus nigra Arn. and Pinus pinaster Ait.) growing in monospecific and mixed stands along the aridity gradient of the Iberian Peninsula show that the productivity of all the studied species, although with different intensity, is modulated by the aridity, which holds with the second hypothesis. The lower the aridity, the higher the productivity, which highlights the importance of taking this variable into account when analyzing productivity and agrees with the first hypothesis Overall, it was found that for each species, the productivity was lower in mixed stands than in monospecific stands (underyielding). This underyielding was, in turn, modulated by aridity. The most evident influence of aridity on mixing effects was found for Ph in the Ph-Pt admixtures and for Pt growing in Pt-Ps mixtures, with both these species presenting overyielding in more arid sites, while in less arid locations the effect changes to underyielding. We only found overyielding for one of the studied species-Pn growing in Pn-Pt mixtures-although this overyielding was also modulated by aridity. These results are consistent with the third hypothesis.
The models obtained allow us to calculate the growth of the species studied and update the volumes of monospecific (and mixtures of two species) pine forests to obtain a national value for a specific year. Through this procedure, it is possible to provide more precise responses to international requirements between cycles of SNFI.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/10/5/430/s1, Table S1: Coefficients resulting from models for monospecific stands without considering mixing effects (Equation (1)), with standard error in parenthesis. Table S2: Comparison between models with and without mixing effect for all pine species. Equation (1) corresponds to the model which does not consider mixtures and Equation (2) considers the mixing effect. Table S3: Variables needed to correct the bias induced by logarithmic transformation (sigma model) and the negative increment of volume (min (VGEsp)). Table S4: Summary of correlations between the principal variables which are used in the models.
Author Contributions: A.A., S.C and M.d.R. conceptualized and designed the methology. A.A. and S.C. developed models and analysis. All authors discussed the results and revised the manuscript.
Funding: This research received no external funding.