Species Mixing Proportion and Aridity Inﬂuence in the Height–Diameter Relationship for Different Species Mixtures in Mediterranean Forests

: Estimating tree height is essential for modelling and managing both pure and mixed forest stands. Although height–diameter (H–D) relationships have been traditionally ﬁtted for pure stands, attention must be paid when analyzing this relationship behavior in stands composed of more than one species. The present context of global change makes also necessary to analyze how this relationship is inﬂuenced by climate conditions. This study tends to cope these gaps, by ﬁtting new H–D models for 13 different Mediterranean species in mixed forest stands under different mixing proportions along an aridity gradient in Spain. Using Spanish National Forest Inventory data, a total of 14 height–diameter equations were initially ﬁtted in order to select the best base models for each pair species-mixture. Then, the best models were expanded including species proportion by area (m i ) and the De Martonne Aridity Index (M). A general trend was found for coniferous species, with taller trees for the same diameter size in pure than in mixed stands, being this trend inverse for broadleaved species. Regarding aridity inﬂuence on H–D relationships, humid conditions seem to beneﬁciate tree height for almost all the analyzed species and species mixtures. These results may have a relevant importance for Mediterranean coppice stands, suggesting that introducing conifers in broadleaves forests could enhance height for coppice species. However, this practice only should be carried out in places with a low probability of drought. Models presented in our study can be used to predict height both in different pure and mixed forests at different spatio-temporal scales to take better sustainable management decisions under future climate change scenarios.


Introduction
Mixed stands have been evidenced as a smart strategy to adapt forest to climate change since they have been shown to provide positive productive and ecosystem service outcomes compared with monocultures [1][2][3][4][5]. Here, diversity plays a significant role in tree growth, environmental adaptation and improvement of meeting demands for goods and services [6][7][8]. Interactions between species may vary with climatic conditions [9] and, hence, climate-growth sensitivity should be explicitly taken into account in the modelling process to consider the impact of climate change [10]. Hence, development of models that consider known species-mixing effects on tree and stand productivity and the potential to improve the resiliency of forestry under expected climate change is essential to properly design the initial composition and subsequent spatial arrangement and management for mixed-species stands [11].
In this context, estimating diameter at breast height and total tree height is fundamental to both developing and applying many growth and yield models [12,13]. The heightdiameter relationship (H-D) allows us to describe stand characteristics and development over time, mean height estimation, stand stability, site index, and growth [14]. Many height-diameter equations have been developed using only diameter at breast height as the independent variable for estimating total height [15][16][17][18][19][20]. However, the height-diameter relationship varies from stand to stand [14], being significantly influenced by other tree and stand characteristics, such as tree health and vigor, site quality and stand development stage, stand density or competition, and species mixture in a stand [12,21,22]. We refer to this type of model as "generalized H-D model" [23]. Thus, the inclusion of additional explanatory variables in the H-D models would make the H-D relationships generalizable over large areas [12,14,[24][25][26]. In Europe, generalized height-diameter models have been fitted and used since the 1930s [27][28][29][30][31][32][33][34], though there are still a limited number of studies focused on tree H-D relationship in mixed forests [22,25,[35][36][37][38]. In this line, recent studies evidence the strong influence of stand structure and species composition on allometric relationships [39][40][41]. However, the magnitude of variations in allometric relationships due to species-mixing is a subject of intense debate, with evidence of adverse and positive effects that depend on multiple factors operating simultaneously [2,39]. 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 [37,42].
On the other hand, climatic conditions are considered one of the main drivers in the biodiversity-productivity relationships showing a wide range of effects across forests globally [5]. However, for instance resistance and relative resilience to drought in mixed stands might depend on depending on species identity, functional diversity in the mixture, site conditions and the developmental stage of the stands [43]. At the same time, in some recent studies, climate variables such as mean annual temperature and precipitation have been found to have a significant effect on tree allometry [44,45]. Hence, H-D models usually designed to be climate-independent will fail to evaluate future scenarios linked to climate change. Previous studies proposed to include of the concept of limiting resources in the tree allometric relationship. Accordingly, [46] added in the modelling approach terms representing abiotic stress, i.e., water availability and [47] included the plant competition symmetry. Tree interactions might be modulated along stress gradients with facilitation being more common in conditions of high abiotic stress (arid sites) relative to more benign abiotic conditions [46]. Therefore, a trend towards higher H-D relationship in mixed stands compared with pure could be expected in the drier and hotter sites, also observed in drought resilience [43].
In the present context of global warming, it is essential to analyze and understand how the H-D relationship is influenced by mixed-species dynamics and environment conditions in order to provide a future management guideline to adapt a Mediterranean forest to climate change. Thus, the inclusion of additional factors on the H-D relationship in mixed forests, such as species proportions and climatic variables, needs to be further considered and analyzed. In this study, we hypothesized that species identity and species composition would define the H-D relationship between in mixed stands compared with pure. We also hypothesized that tree height in mixtures could be benefited under drought conditions compared with pure stands. In this context, this paper focused on analyzing how the identities of competitors and aridity influence the height-diameter relationship in Mediterranean mixed forests. To achieve that, new generalized H-D models that allow us to estimate tree total height in mixed forests along a gradient of different mixtures and species proportions were fit along an aridity gradient in Spain.

Data
In this study, data from the second edition of the Spanish National Forest Inventory (SNFI) were used. The SNFI is a systematic sample where plots are located at the nodes of a kilometer square grid and remeasured every ten years. SNFI plots are composed of four concentric sample circles of 5 m, 10 m, 15 m and 25 m radius, where the diameter and height of all trees over 7.5 cm, 12.5 cm, 22.5 cm and 42.5 cm breast height diameter are recorded, respectively [48]. Stand-level variables were calculated using tree measures and expansion factors. Annual mean temperature (T) and annual total precipitation (P) were derived from Worldclim 2 [49] for the selected SNFI plots (Supplementary Table S1). With the aim of analyzing how species composition influences H-D relationship in mixed stands, a broad range of mixing conditions were considered in our study, selecting both pure and mixed plots (Supplementary Figure S1). We focused on 13 main tree species: Aleppo pine (Pinus halepensis Mill.), Black pine (Pinus nigra Arnold), Maritime pine (Pinus pinaster Ait.), Stone pine (Pinus pinea L.), Scots pine (Pinus sylvestris L.), Pinus mugo (Pinus uncinata Turra), Beech (Fagus sylvatica L.), Portuguese oak (Quercus faginea Lam.), Holm oak (Quercus ilex L.), Sessile oak (Quercus petraea Matt. Liebl.), Pyrenean oak (Quercus pyrenaica Willd.), Pedunculate oak (Quercus robur L.) and Cork oak (Quercus suber L.). First, SNFI mixed plots (composed by two main species) were selected based on the following criteria: (1) the main species are included in the pool of studied species; (2) the proportion of both species (based on basal area) being higher than 90% and (3) the proportion of each main species being higher than 15%. SNFI plots were considered as pure if the proportion of basal area of the main species was 90% or higher. Only pure stands close to mixed stands were taking into account for creating the study dataset in order to guarantee similar site conditions. To achieve that, a buffer of 3 km of diameter was created around the dissolved polygon of mixed plots and pure plots, and final plot selection was made inside the intersection surface. All spatial processes were carried out in QGis software [50]. We discarded all species compositions with less than 50 mixed plots in the dataset to keep similar site conditions, which results in 29 representative species compositions being analyzed.

Base H-D Model Selection
Model fitting was carried out in two steps. In a first step, basic height-diameter (H-D) models were fitted for each species and species compositions studied. A total of 14 expressions of H-D equations most used from the forestry literature (Table 1) were considered. These base candidate models were fitted for each species and mixture composition using the optimize.curve_fit function from the scikit-learn [51] library in Python. To validate models and check overfitting, a cross validation bootstrap process was carried out. Data were split into train and test sets in an 80/20 proportion for each pair of mixture composition. The train set was used to obtain parameter estimates, since the test set was used for model evaluation performance. Confidence intervals for each parameter were obtained by using a nonparametric bootstrap procedure, as described by Robinson and Froese [52]. The number of bootstrap replicates was set to 1000. Following this bootstrap procedure, root mean squared error (RMSE) and Akaike Information criterion (AIC) were obtained for each model. The best H-D base model was selected for each species in a mixture composition based on parameter significance and model performance (lower AIC and RMSE). Then, selected H-D base models were expanded to allow the inclusion of competition and climate variables in the H-D relationships.

Extension of the Base H-D Models
In a second step, the previously selected H-D base models were expanded including several independent variables representing stand development, competition, species proportion and climatic conditions. Since we considered these variables to mostly affect asymptotic parameter (β o ) of the base H-D models, this parameter was further expanded both as a linear (Equation (1)) and multiplicative (Equation (2)) form as follows: where β o is asymptotic parameter of base H-D models presented in Table 1; BAL is the basal area of trees larger than a subject tree, in m 2 ·ha −1 ; Dq is the quadratic mean diameter, in cm; m i is the species proportion by area of species i, described below; M is the De Martonne aridity index, in mm· • C −1 ; and Ho is the dominant height, in m. Note that Dq and Ho were tested both at plot and species (Dq i , Ho i ) level since preliminary analysis showed a high correlation between these variables at species level with height estimation, especially for Quercus species.
Then, the best H-D model was selected for each species and mixture, considering variable significance and model parsimony using AIC corrected for small sample size (AIC c ). After that, we fitted non-linear mixed effects models to consider hierarchical and spatial data structure of forest inventory data. Numerous studies have applied the mixedeffects models to describe H-D relationships and have alleviated this lack of independence of error terms, improving the model fitting and prediction accuracy [53][54][55][56]. To consider random effects in the expanded H-D models, the following procedure was carried out: Firstly, we tested plot as random variable based on a similar study [37]. This variable was included iteratively into the different parameters of the expanded H-D models to determine the best model for each pair of species-species composition to allow their comparison, models were fitted using maximum likelihood (ML). Secondly, among the pool of fitted models, we identified the best expanded mixed-effects H-D model for each pair species-mixture composition considering model parsimony using an informationtheoretic approach. The model with the lowest AIC c and greater Akaike weight (W i ) was considered the best and most parsimonious model for the observed data relative to the set of alternative models. Finally, the best-fitting expanded H-D functions were re-fitted using the unbiased restricted maximum likelihood method (REML) with R-package 'nlme' [53]. To ensure that our modelling strategy has accounted for heteroscedasticity, models were assessed using residuals plots. We assessed the contribution of forest stand and climate variables by looking at the significance of their respective parameters.

Mixing Proportion Influence on the H-D Relationship in Mixed Forest Stands
We used species proportion by area (m i ) to characterize stand composition and to identify the significant effects of species interactions on H-D relationships in mixed stands. The species-specific growing space occupied is relevant for calculation of the mixing proportions, stand density and quantification of mixing effects on growth [57,58]. Thus, species mixing proportions might be calculated to avoid bias in the quantification of the net total mixing effect, as well as in the relative importance of under-or overyielding by species, due to differences in the potential densities of the species [59]. Competition equivalence coefficient (CEC) compares species-specific growing space requirements of a species with their value in mixed stands, and it is calculated as the ratio of potential carrying capacity of both species in pure stands (by means of maximum stand density index [60]-SDI max ). Since recent studies have shown evidence that SDI max varies with climate [61,62], in this study we obtained climate-dependent CECs (Supplementary Table S2) to calculate species mixing proportions for each species in mixed stands using climate-dependent Maximum Size-Density Relationship (MSDR) models presented by Rodriguez-de-Prado et al. [63]. Table 1. Base tree Height-Diameter equations used in this study. h = total tree height; d: diameter at breast height; Ho = dominant height of the plot; do = dominant diameter of the plot, β 0 , β 1 , and β 2 : model parameters.

Model
Reference Formula Therefore, species proportion by area was calculated for each species based on SDI and equivalence coefficients following similar studies [57,[77][78][79]: where m i is the species proportion by area of species i; SDI i is the stand density index of species i; SDI j is the stand density index of species j; CEC j,i is the competition equivalence coefficient for species j, taking as reference species i; SDI max,i is the climate-dependent maximum stand density index of species i; and SDI max,j is the climate-dependent maximum stand density index of species j.
To compare the magnitude of the mixture effect on tree height, the M (ratio) was calculated, following a similar approach to that of [41]. For each of the studied species and mixtures, tree height for a given diameter variation according to the species competition environment was estimated for pure stands and mixed stands using the fitted models for a reference tree of 30 cm diameter, and the other independent variables were fixed to the mean (Supplementary Table S1). M (ratio) was calculated as follows: where H (mix) is the estimated mean total height of species i in mixture j, at a species mixing proportion m i , and H (pure) is the estimated mean total height of species i in mixture j at m i = 1. Thus, we quantified the relative change on tree height for a tree of 30 cm diameter in response to inter-specific competition along a gradient of different species mixing proportions (m i from 0.2 to 0.8). Therefore, when M (ratio) > 1, the estimated total tree height would be higher in mixed than in pure stands, and lower tree height is expected in mixed stands when M (ratio) < 1.

Climatic Influence on H-D Relationships in Mixed Forests Stands
The De Martonne Aridity Index (M) was included in the models as a measure of aridity (M = P/(T + 10); mm • C −1 ) [80]. Including this index into H-D models allowed us to estimate total tree height along an aridity gradient. Accordantly with M (ratio) , C (ratio) was here calculated to analyze potential differences of total tree height for a given tree diameter (30 cm) between arid and humid environments. All independent variables (excepting M, and m i = 0.5) were fixed to the mean values for each species and mixture (Supplementary Table S1). C (ratio) was proposed and calculated in this study as follows: where H (arid) is the total tree height estimated using fitted models for species i in mixture j along the aridity gradient, from smallest to the highest values of the De Martonne Index. H (humid) is the estimated total tree height predicted at the highest De Martonne Index found. C (ratio) depicts the magnitude of relative tree height for a given diameter higher (C (ratio) > 1) or lower (C (ratio) < 1) in more arid than humid places, respectively.

New H-D Models for Mixed Forests in Spain
Coefficients for the final best generalized mixed-effects H-D models for each species and mixture are presented in Table 2. Results indicated that M1 [64] was the most frequently chosen base H-D model among conifers, whereas M12 [74] was the most chosen among broadleaves. An exception was found for Pinus sylvestris and Pinus uncinata, where M3 [66] and M5 [68] showed better fits for H-D relationships, respectively. All the parameters of the expanded H-D models obtained by bootstrap procedure were significant (p-value < 0.05), finding some differences regarding the independent variables entering the models for coniferous and broadleaves. For the coniferous species, competition (BAL) and species proportion by area (m i ) had a significant effect explaining height in the models. For these species, the competition at plot level (by means of Dq) was included in the models more times when mixed with conifers, but Dq i (at species level) was included more often when mixed with broadleaves. Indeed, for all the Quercus species analyzed, Dq i and Ho i were included into the H-D expanded models as independent variables. Note that, for almost all the species and species compositions of broadleaved species, the De Martonne Index (M) was significant with a positive sign, meaning a height increment by humid conditions. For all pairs of species mixtures, expanded H-D models significantly improved the goodness of fit, in terms of AIC and RMSE, compared with the base H-D models. Among the analyzed species, Pinus pinea and Quercus ilex selected models showed the smallest RMSE with values ranging 0.80-1.06 m. On the contrary, the highest RMSE values were found for Pinus sylvestris, Pinus uncinata, Quercus petraea, Quercus pyrenaica and Quercus robur with values close to 1.5 m.   Figure 1 shows the height variation at the same diameter size under different species mixing proportions in mixed forest stands for the studied species. A general trend was found for Pinus species (Figure 1a-f), with taller trees in pure than in mixed stands (M (ratio) < 1). However, an exception was found in for Pinus sylvestris in Pinus sylvestris-Pinus pinaster mixture (Figure 1e), with slightly higher trees at the same diameter size in the mixture than pure stands of this species. Pinus uncinata height was nearly insensitive to proportional area change and mixed and pure forest allometry (Figure 1f). On the other hand, an inverse trend was found for broadleaved species, where trees were comparatively taller (~2-8%) in mixed stands than in pure ones (Figure 1g-m) at the same diameter size. Among them, Fagus sylvatica, Quercus faginea and Quercus robur height was comparatively stable under different species mixing, finding similar estimated tree heights both in pure and mixed forests (Figure 1g,h,l). Among the conifers, the highest differences in height between mixed and pure stands were observed for Pinus nigra (Figure 1b) when mixed with Quercus faginea and Quercus ilex (~10% reduction in height in mixtures). Regarding broadleaved species, Quercus ilex showed a high difference between mixed and pure stands when mixing with Pinus halepensis (Figure 1i), with taller trees at the same diameter in mixed stands. This trend was also found for the Quercus suber and Pinus pinaster mixture ( Figure 1m); Quercus petraea also presented significant differences in height with respect to pure stands when mixed with Fagus sylvatica (Figure 1j).

Aridity Influence on H-D Relationships
The variation of height at the same diameter size under different aridity conditions, by terms of the De Martonne Aridity Index, is presented in Figure 2. According to our results, taller trees could be found in more humid conditions for almost all the analyzed species and species mixtures. Based on estimated height along an aridity gradient, results also indicated that conifers studied here are less sensitive to changes in aridity than broadleaved species (Figure 2). Among the studied species and species compositions, aridity influence was nearly insignificant for Pinus nigra, Pinus pinaster and Quercus suber (Figure 2b,c,m). Two different trends were found regarding climatic influence on height for conifers. On one hand, Mediterranean conifers such as Pinus halepensis (Figure 2a) and Pinus pinea (Figure 2d) are approximately 5% higher in more arid than humid conditions (C (ratio) < 1) when mixed with Quercus ilex, mainly. Oppositely, estimated heights for species living at higher altitudes such as Pinus sylvestris (Figure 2e) and Pinus uncinata (Figure 2f) presented higher values at more humid environments (C (ratio) > 1), although these differences were nearly insignificant (~1-2%). Broadleaved species (Figure 2g-m) seemed to be highly influenced by aridity according to their height estimations, excepting Quercus suber. Among them, the highest differences in height between arid and humid conditions were found for Quercus faginea in the mixture with Pinus sylvestris (Figure 2h) and Quercus robur in mixture with Fagus sylvatica (Figure 2l), with values close to 15% (height difference in humid compared drought sites). Figure 3 shows total tree height estimations (y-axis) along the diameter distribution of each species in each mixture (x-axis) under different values of dominant height. Here, we can see that trends in aridity and height found for each species and mixture are constant under different tree sizes and stand developmental stages. In the case of Pinus halepensis, taller trees were found in more arid than in less arid places for all sizes and developmental stages (Figure 3a) when mixed with Quercus ilex. On the other hand, height estimations for Quercus robur when mixed with Fagus sylvatica indicated that, under the three simulated stand developmental stages, taller trees were found in more humid places than in more arid places (Figure 3b).  y-axis). Letters a-m indicate the species name in order to make reference to them along the text. Note: M (ratio) represents whether height in mixed forests is higher (M (ratio) > 1) or lower (M (ratio) < 1) than pure stands. All independent variables (excepting m i ) have been fixed to the mean of each species.  Letters a-m indicate the species name in order to make reference to them along the text. Note: C (ratio) ratio (y-axis) represents whether height in mixed forests is higher (C (ratio) > 1) or lower (C (ratio) < 1) in more arid than humid places. All independent variables (excepting M and m i = 0. 5) have been fixed to the mean values for each species.

Figure 2.
Variation of height along the aridity gradient (x-axis), with Arid being the lowest M value and Humid the highest M found in the distribution of each species in each mixture composition. Letters a-m indicate the species name in order to make reference to them along the text. Note: Cratio ratio (y-axis) represents whether height in mixed forests is higher (Cratio > 1) or lower (Cratio < 1) in more arid than humid places. All independent variables (excepting M and mi = 0. 5) have been fixed to the mean values for each species.

Discussion
This study presents novel non-linear height-diameter models for predicting tree height in Mediterranean mixed forests under different species mixing proportions and aridity conditions. Here, we evidenced that stand composition, target species, identity of competitors and aridity conditions modulate H-D relationship for the species and mixtures analyzed, resulting in positive, neutral or negative effects on H-D relationships. Compared with existing models to estimate tree height in Spain [81], our work contributes towards H-D models accounting for aridity conditions and species composition in the estimates in addition to tree competition and plot-level factors. Tree height estimates accurately represent species mixing proportions and aridity conditions will improve derived upscaling volume estimations from tree to stand level. Hence, our results could be important to support management and policy decision for Mediterranean forests under the context of climate change, especially for mixed stands.

Discussion
This study presents novel non-linear height-diameter models for predicting tree height in Mediterranean mixed forests under different species mixing proportions and aridity conditions. Here, we evidenced that stand composition, target species, identity of competitors and aridity conditions modulate H-D relationship for the species and mixtures analyzed, resulting in positive, neutral or negative effects on H-D relationships. Compared with existing models to estimate tree height in Spain [81], our work contributes towards H-D models accounting for aridity conditions and species composition in the estimates in addition to tree competition and plot-level factors. Tree height estimates accurately represent species mixing proportions and aridity conditions will improve derived upscaling volume estimations from tree to stand level. Hence, our results could be important to support management and policy decision for Mediterranean forests under the context of climate change, especially for mixed stands.

Total Tree Height Response to Species Mixing Proportions
The information on tree heights obtained from H-D models is essential in forest inventories for computing tree volumes and evaluate tree slenderness which could be related to stand stability [82]. Our results showed that tree volume seemed to be higher in pure stands for conifer, but with lower stability (high slenderness), increasing the risk of damage in windstorms. On the contrary, slenderness and tree volume increased in mixtures for broadleaves species, which could be key to stimulating tree development for Mediterranean coppice stands. Taking into account species mixing proportions in mixed forest studies is key to understanding differences in tree and stand variables of a given species [79]. In this line, our results suggested that conifers studied here could be taller at the same diameter size in pure than in mixed forest stands (Figure 1a-f). We hypothesized that it may be caused by species-specific traits, such as pines are pioneer and shade intolerant species; thus, they are very sensitive to growing under broadleaves cover. This depends primary on species identity, i.e., shade tolerance is caused for plant ontogeny and influenced by numerous biotic and abiotic factors [83]. Shade tolerance represents how species intercept light and how efficiently it is converted into photosynthates; hence, we would expect that mixtures of different shade tolerance could affect tree light competition, growth and H-D relationship. Accordingly, [84] showed that species-specific shade tolerance may affect tree allometry and be age-mediated.
Among the conifers studied here, Pinus nigra showed the biggest differences in height between mixed and pure stands (Figure 1b), especially when mixed with oaks. It is known that although this is a high-shade-tolerant species compared with other Mediterranean pines, water limitation could drastically limit its distribution in Mediterranean plant communities [85,86]. An exception among conifers studied here was found in Pinus sylvestris mixed with Pinus pinaster, with higher trees at the same diameter size found in mixed than in pure stands, although the M (ratio) keeps close to 1 along all the mixing proportions gradient (Figure 1e). A similar trend was found by Riofrio et al. [37] for this species and mixture composition from SNFI data. Another related explanation to these patterns could be due to differences in species growth rates. When the mixture occurs between two species with different growth rates, the estimated height for the pure stand may be higher at the same stage of development [87]. This seems logical because all trees in a regular pure stand grow very similarly and in a regular way, maintaining the upper average and competing for light availability. This is much more evident for pine in conifer-broadleaved mixtures, where the latter do not grow as much in height and the crown structure is not pyramidal [36], increasing inter-specific competition for pine species. Competition may also explain our results found for Quercus petraea when mixed with Fagus sylvatica (Figure 1j), with great differences in height at the same diameter size between mixed and pure stands for oak species. In mixture, we hypothesized that the growth rate of the light-demanding oak is supposed to be higher than the shade tolerant [3]. Thus, an increase of beech proportion in mixture could reduce intra-specific competition for oak and enhance height for this species. However, this contradicts Ligot et al. [88] observations in forests of Central Europe, where oak trees were systematically outcompeted by beech in the mixture, preventing also oak regeneration. Focusing on broadleaved species, and contrary to conifers, a common trend was found where higher trees at the same diameter size were found in mixed stands (Figure 1g-m). Among other reasons, this inverse trend could be due to stand stratification promoted by silvicultural history. Silviculture can accelerate the growth of different species in mixtures, affecting stand composition, structure and dynamics [89]. In some stands, one of the species in the mixture may historically experience repeated fires, grazing or silvicultural treatments (cutting) to favor the other species. This can be the case of the Quercus-Pinus mixtures analyzed, where H (mix) /H (pure) ratios higher than 1 for broadleaved in mixed stands may suggest a positive effect of the mixture. Note that the proposed M (ratio) index relied on species plasticity in tree allometry in response to competition, which, according to Pretzsch (2014) [2], can be defined as the capacity of the species to change tree allometry under different competitive environments. As broadleaf coppice forests are wide extended in Mediterranean areas, we think that the mixture with conifers may reduce the intra-specific competition and thus may facilitate the apparition of taller trees. In this sense, our results may have a relevant importance for Mediterranean coppice stands, suggesting that introducing or combing conifers in broadleaved forests, could enhance the height-diameter relationship, obtaining taller trees at the same diameter for coppice species. This practice would also increase the ecosystem services and productivity; thus, it has been commonly used in reforestation [90]. Differences in trends between conifers and broadleaved species studied here could be also explained based on initial SNFI selected plots, which could have different stand structures and development stage of the species, as shown in Supplementary Figures S3 and S4. It seems that in conifer mixtures, the size of the trees seems to be different in strict Mediterranean species (Pinus halepensis, Pinus pinaster, Pinus pinea) than in montane species (Pinus sylvestris, Pinus nigra, Pinus uncinata), probably as a consequence of the character of greater shade intolerance. In addition, as the stand distribution for a species moves away from its ideal ecological niche, competition with other more adapted and opportunistic species appears. If the second or third species is equally or better adapted, it leads to the formation of mixed stands. This can be the reason why Quercus species may perform better in mixtures (mainly broadleaved-conifer as Quercus ilex-Pinus halepensis or Quercus suber-Pinus pinaster) than pines, since they naturally form this kind of mixed forests. Differences in water uptake depth appears to be stimulated in mixed forests because Quercus ilex surrounded by Pinus halepensis explores deeper water sources than in monospecific formations [91]. We observed that shade-tolerant species as Fagus sylvatica, Quercus robur and Pinus uncinata, were no affected by mixture, regardless species composition. This may be due to shade tolerant species having high plasticity for certain traits, particularly for morphological features which optimize the light capture [83]; therefore, a high canopy cover promoted by mixed stand seems to be practically irrelevant.

Total Tree Height Response to Aridity
We observed a significant climatic influence on the H-D relationship for most of the species under study, in line with similar research in the topic [44]. Including climatic related variables in the H-D models could be the key to design forest prescriptions under different present and future environmental conditions [92]. A basic assumption in tree allometry studies is that, as long as site conditions are homogenous, trees of a given species, at a given location, with the same diameter at breast height, would have the same height [78]. However, climatic conditions change over locations and time, and they could contribute towards species-specific dynamics [93,94]. In this context, it is well known that height-diameter relationships depend heavily on local environmental conditions [95]. In the present scenario of climate change, the fluctuating climatic conditions may have a positive or negative influence on the H-D relationship [44]. We observed that for most of the mixtures analyzed here, total tree height along a diameter range was reduced under drought conditions ( Figure 2). However, two exceptions opposed to this trend were observed: Pinus halepensis-Quercus ilex (Figure 2a) and Pinus pinea-Pinus pinaster (Figure 2d). Although in the first case it may be due to a reduction in intra-specific competition for Pinus halepensis during drought conditions, in the second case, it could be due to a facilitation process, i.e., a complementarity in the use of resources [96]. Due to its stronger plastic character and its potential adaptability, Pinus halepensis has been demonstrated to be the most suitable species in terms of tree growth in arid sites of Spain [97]. A more efficient use of below-ground water resources due to different root stratification could be also speculated for these mixtures, which allows for a better exploration of the soil profile [91,98]. Therefore, these species could use more resources to strengthen the root system than the aerial system as aridity increases.
Potential changes in height due to changes in aridity may be higher in Mediterranean broadleaved than conifers studied here, suggesting that the allometry of oak species in Spain may be highly affected by future climate change. We hypothesized that the differences between species on water use efficiency may underlay this process. Pinus species can maintain a relatively stable leaf water potential by strict stomatal control during drought events (isohydric strategy). On the other hand, broadleaves do not have a discernible threshold of minimum water potential response, i.e., a light stomatal control under drought conditions (anisohydric behavior) [96][97][98][99]. Consequently, higher competition for water use resources in arid sites is expected for broadleaves when mixed with conifers [100]. Our results are in line with this hypothesis, as shown especially for Quercus faginea when mixed with Pinus sylvestris (Figure 2h). As previously observed in other oak-pine mixture compositions [101], the impact of soil water deficit on species may be worsened in mixtures, but a mixture of species with different growth sensitivities to the seasonality of the drought periods might help to buffer these effects.

Conclusions
Data from the Spanish National Forest Inventory (SNFI) were used to estimate the H-D relationship, covering thirteen species in twenty-nine different species compositions in Spain. The modified generalized mixed model was more adequate for modeling the height-diameter relationship of the tree species groups over the base H-D models. We found a significant contribution of climate and stand species composition variables in modifying the H-D relationship. Regarding species mixing proportions, taller trees at the same diameter size were found more commonly in monospecific stands for conifers and mixed stands for the broadleaved species studied here. Based on aridity conditions, in general, taller trees at the same diameter are supposed to be found in more humid conditions, with exceptions found for Pinus pinea and Pinus halepensis species in speciesspecific mixtures. The broadleaved and conifer mixed stands showed different patterns. Broadleaves species are more sensitive to prone drought sites. The use of forest inventory data allowed us to take into account a wide range of different species and conditions (climatic, soils, competition, species composition or mixtures). H-D models presented in this study permit us to save time and effort by predicting tree heights instead of performing direct measurements. The use of new programming techniques such as machine learning and deep learning applied to more precise data (LiDAR and satellite imaginery) could help us obtaining more precise height estimations in the future within the scope of further research. H-D models presented in our study can be an important tool to estimate height at different spatio-temporal scales in order to take better sustainable management decisions under future climate change scenarios.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/f13010119/s1: Table S1: Stats (mean and range) of the main tree and plot variables characterizing selected SNFI plots used in our study; Table S2: Competition Equivalence Coefficients for the most representative species compositions in Spain; Figure S1: Location of the Spanish National Forest Inventory plots selected to fit the H-D models for the 29 different mixtures analyzed in this study; Figure   Acknowledgments: The authors would like to thank Susana Miguel Barreiro for revising the complete text. In addition, authors would also thank to MSCA-RISE "3D-NEONET: Drug Discovery and Delivery NEtwork for ONcology and Eye Therapeutics" project for the opportunity of working close to the National University of Galway, School of Computer Science in the development and use of new programming techniques applied to sustainable forest management. Finally, the authors would like to thank to the Ministry for Ecological Transition and the WorldClim team for sharing and providing the data used in this study.