National-scale assessment of forest site productivity in Spain

Contents lists available at ScienceDirect Forest Ecology and Management journal homepage: www.elsevier.com/locate/foreco National-scale assessment of forest site productivity in Spain Daniel Moreno-Fernández a,b,⁎ , Juan Gabriel Álvarez-González c , Roque Rodríguez-Soalleiro d , María Pasalodos-Tato a , Isabel Cañellas a , Fernando Montes a , Emilio Díaz-Varela e , Mariola Sánchez-González a , Felipe Crecente-Campo f , Pedro Álvarez-Álvarez g , Marcos Barrio-Anta g , César Pérez-Cruzado c a INIA-CIFOR, Ctra. A Coruña km 7.5, E-28040 Madrid, Spain b MONTES (School of Forest Engineering and Natural Resources), Universidad Politécnica de Madrid, 28040 Madrid, Spain c Department of Agroforestry Engineering, Unit of Sustainable Forest Management (UXFS), School of Engineering, University of Santiago de Compostela E-27002 Lugo, Spain d Department of Forest Production and Engineering Projects, Unit of Sustainable Forest Management (UXFS), School of Engineering, University of Santiago de Compostela, E-27002 Lugo, Spain e Department of Plant Production and Project Engineering, ECOAGRASOC Research Group, School of Engineering, University of Santiago de Compostela, E-27002 Lugo, Spain f CERNA, Illas Cies 52, E-27003 Lugo Spain g GIS-Forest Research Group, Department of Organisms and Systems Biology, Polytechnic School of Mieres, E-33600 Mieres, Asturias, Spain ARTICLE INFO Keywords: Uneven-aged stands Height-diameter curves Sustainable forest management International forest policy Guide curves ABSTRACT Sustainable production of wood is one of the main services provided by forest systems. Site productivity in the case of forests is often evaluated through the site quality. However, most of the works addressing the site quality have been done at local or regional scale. In this work, we aim to develop site quality models for five dominant species in Spanish forests (Fagus sylvatica,Pinus pinaster atlantica,Quercus pyrenaica,Pinus nigra,Pinus sylvestris) and create site quality maps at a national-scale from these models. First, we develop site quality models using site form (height-diameter relationship) as the reference index and the Spanish National Forest Inventory as dataset. Then, we fit spatial additive models entering physiographic and climatic variables in order to predict the site quality over the whole country. Additionally, we plot site form maps for the five species in order to describe spatial pattern in site quality at a national scale. Altitude and aspect appeared to be fundamental variables in the assessment of site quality. The accuracy of the spatial additive models ranged from 38.2% to 47.9%. The correspondence between the predicted and observed maps of site qualities is clear. Our results provide a tool which could be used by forest managers in land use planning as well as in forest policy decision-making at a national scale. We suggest that this method could be used in other countries and that the maps could be expanded to the European scale to assessing the way in which site quality varies across Europe always considering that the relationships between forest productivity and environmental variables could vary among biogeoclimatic zones. 1. Introduction Site quality has long been used in forestry as a proxy of site productivity. This variable is not only of interest for predicting growth and yield of forest stands (Álvarez-González et al., 2005; Clutter et al., 1983; Diéguez-Aranda et al., 2005) but also for studies on ecological diversity (Franklin et al., 1989), forest structure (Larson et al., 2008) and forest disturbances (Wei et al., 2003) among others. Site quality has traditionally been expressed as the relationship between dominant height and age (site index) for even aged forests. However, this definition is difficult to apply in forests with an uneven distribution of ages or in mixed-species forests, where the height-diameter relationship (site form) has proved to be a good measure of site productivity (Huang and Titus, 1993; McLintock and Bickford, 1957; Stout and Shumway, 1982; Vanclay, 1994). Site form can be also useful to determine the site quality of even-aged stands where the age is unknown, as occurs in the National Forest Inventories (NFIs) of some countries (Tomppo et al., 2010). https://doi.org/10.1016/j.foreco.2018.03.016 Received 11 January 2018; Received in revised form 6 March 2018; Accepted 9 March 2018 Abbreviations: AIC, Akaike’s Information Criterion; B-R, Bertalanffy-Richards model; dbh, diameter at breast height; EF, expansion factor; EFT, European Forest Type; H-II, Hossfeld II model; NFI, National Forest Inventory ⁎ Corresponding author at: INIA-CIFOR, Ctra. A Coruña km 7.5, E-28040 Madrid, Spain. E-mail address: [email protected] (D. Moreno-Fernández). Forest Ecology and Management 417 (2018) 197–207 Available online 16 March 2018 0378-1127/ © 2018 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/BY/4.0/). T Site form, as well as site index, can be estimated via direct and indirect methods. Direct methods are based on the relationship between the height and diameter, in the case of site form, and between height and age, in the case of site index. Meanwhile indirect methods relate physiographic, climatic, edaphic variables and understory to site quality (Bravo and Montero, 2001; Pacheco Marques, 1991). Alternative base models have been employed to develop diameter-height relationships. Meyer’s mathematical expression have been used for fitting site form curves in Picea rubens Sarg. (McLintock and Bickford, 1957) and in six hardwood species (Stout and Shumway, 1982) while Huang and Titus (1993) selected the Bertalanffy-Richards model as the base height-diameter function. Different statistical methodologies have been used to determine the influence of environmental variables on the site quality, such as linear models (Chen et al., 2002; Pacheco Marques, 1991; Seynave et al., 2005), discriminant analysis (Bravo-Oviedo and Montero, 2005; Bueis et al., 2016) or regression trees (Álvarez-Álvarez et al., 2011). The group of variables selected for site quality models and those variables found as significantly associated to site quality have varied among species and study sites. Bravo-Oviedo et al. (2011) stated that variations in site quality of Mediterranean Pinus pinaster Ait stands are mainly explained by climate variables. Other authors, however, only considered edaphic and physiographic variables in site quality models (Bravo et al., 2011; Bueis et al., 2016; Seynave et al., 2005). Aspect, latitude, content of nutrients, soil moisture, soil texture, pH or soil depth appeared to be related to site quality (Bravo and Montero, 2001; Bueis et al., 2016; Seynave et al., 2005). Nevertheless, the effects of a variable on the site quality of a given species may change among regions (Chen et al., 2002). Additionally, some of these relationships might not be linear as in the case of the site quality-altitude (Seynave et al., 2005). In these cases, linear models cannot identify complex relationships among variables. Therefore, approaches, such as additive models, describing nonlinear and complex relationships between the site quality and the predictors can be very useful (Hastie and Tibshirani, 1989; Wood, 2006). Most of the studies assessing the quality site have been carried out at local or regional scales (Bueis et al., 2016; Seynave et al., 2005) whereas the knowledge of forest productivity at larger geographic scales is scarce (Chen et al., 2002). In this regard, some authors affirm that the representativeness at a coarser scale of forest productivity from restricted area data is questionable (Charru et al., 2010). In this sense, NFIs provide the broadest source of knowledge on the status of the forests at national level in many countries (Barbati et al., 2014). In fact, NFIs have been used to estimate production (Charru et al., 2010), aboveground biomass (Avitabile and Camia, 2018), carbon storage (Woodall et al., 2008) and monitoring species distribution (Hernández et al., 2014; Moreno-Fernández et al., 2016) or assessment of biodiversity (Andersson and Östlund, 2004) and abiotic damages (Jalkanen and Mattila, 2000). Additionally, the NFIs have been used as datasets in scientific studies to model recruitment (Lexerød, 2005), tree biomass (Ruiz-Peinado et al., 2011), tree mortality (Ruiz-Benito et al., 2013), deadwood volume (Crecente-Campo et al., 2016) as well as site quality (Adame et al., 2006). International institutions and processes, however, favor using Table 1 Number of plots of the National Forest Inventory (N), number of plots with basal area of the species in parenthesis larger than 90% (N G ≥ 90%), mean basal area (G in m 2 ha −1 ), mean dominant height (H 0 in m) and mean dominant diameter (D 0 in cm) in the selected European Forest Types. Minimum and maximum values in brackets. European Forest Type N N G ≥ 90% G H 0 D 0 2.7 Atlantic maritime pine forest (Pinus pinaster atlantica) 3410 1568 18.0 (0.4–85.0) 15.2 (2.0–29.3) 30.5 (7.5–71.9) 7.1 South western European mountainous beech forest (Fagus sylvatica) 4449 1915 25.5 (0.4–64.2) 18.8 (3.7–36.5) 40.9 (7.5–158.8) 8.3 Pyrenean oak forest (Quercus pyrenaica) 5528 2889 10.9 (0.4–72.5) 9.9 (0.4–23.3) 24.2 (7.5–134.4) 10.2 Mediterranean and Anatolian black pine forest (Pinus nigra) 8352 3212 15.1 (0.4–80.4) 10.4 (2.6–33.3) 26.5 (7.5–100.3) 10.4 Mediterranean and Anatolian Scots pine forest (Pinus sylvestris) 10,919 5171 22.7 (0.4–83.8) 12.1 (2.0–30.5) 29.9 (7.5–64.8) Table 2 Base model, site specific parameter and site specific parameter for each site form family of curves. Species Base model Parameter Equation Estimate Pr(> |t|) P. pinaster H-II a b = 0.13035 < 2 10 −16 = + + SFI 1.3 Dref D Hb Doref D 0 2 0 01.3 (0) 2 F. sylvatica H-II a = 2.5473 < 2 10 −16 = + + SFI 1.3 Dref aD Ha Dref D 0 2 0 01.3 0 0 2 Q. pyrenaica H-II a = 2.0539 < 2 10 −16 = + + SFI 1.3 Dref aD Ha Dref D 0 2 0 01.3 0 0 2 P. nigra B-R b b = 0.05405 < 2 10 −16 = +SFI H1.3 ( 1.3) ebD ref ebD c 010 10 c = 2.2316 < 2 10 −16 P. sylvestris H-II b = 0.1386 < 2 10 −16 = + + SFI 1.3 Dref D Hb D ref D 0 2 0 01.3 (00) 2 a H-II: Hossfeld-II. b B-R: Bertalanffy-Richard. D. Moreno-Fernández et al. Forest Ecology and Management 417 (2018) 197–207 198 international rather than national criteria (FAO, 2015; Gabler et al., 2012) to monitor sustainable forest management practices and forest state at international scales (Larsson, 2001). In this context, Barbati et al. (2014) and Barbati et al. (2007a) have proposed a classification to categorize the forest of the Pan-European region into 14 European Forest Types (EFTs) categories (category level) and the categories into 78 EFTs (type level) according to ecological sound units. The EFTs classify forest features by aggregating and averaging data from NFIs into ecologically homogeneous strata of European relevance, i.e. the 14 EFTs categories (Barbati et al., 2014). Therefore, reporting forest productivity through an index independent of the forest structure and following EFTs classification would contribute to standardize assessment resources in Europe or even abroad. In this work, we aim to develop site quality models at a nationalscale using site forms for five species, and display the results spatially. In the first step, we fit site form models (diameter-height curves) for the Fig. 1. Site form family of curves for the European Forest Types studied. Table 3 p-values of the model components and percentage of deviance explained by the spatial additive model for each European Forest Type. Species Intercept f(X i ,Y i )f(altitude)f(aspect)f(slope) Deviance (%) P. pinaster < 0.0001 < 0.0001 < 0.0001 n.s. a n.s. 39.6 F. sylvatica < 0.0001 < 0.0001 < 0.0001 < 0.0007 0.0400 37.3 Q. pyrenaica < 0.0001 < 0.0001 < 0.0001 0.0192 n.s. 32.5 P. nigra < 0.0001 < 0.0001 < 0.0001 < 0.0001 n.s. 48.6 P. sylvestris < 0.0001 < 0.0001 < 0.0001 < 0.0001 n.s. 42.0 a n.s. = non-significant. D. Moreno-Fernández et al. Forest Ecology and Management 417 (2018) 197–207 199 five species using NFI as dataset. We then model the site form using spatial additive models including environmental variables as predictors. Finally, we used these spatial additive models to create site quality maps by species at a national-scale. 2. Material and methods 2.1. Species studied and National Forest Inventory dataset We used the EFT classification to select the target species. We selected five species widely distributed in Spain and which form genuine, monodominant forest according to the EFT classification (Barbati et al., 2014, 2007b): Pinus pinaster Ait. ssp atlantica (P. pinaster, hereafter), Fagus sylvatica L., Quercus pyrenaica Willd., Pinus nigra Arn. and Pinus sylvestris L. (Table 1). Hereafter, each EFT is named after the dominant species present. All of the species are found mainly in montane areas with the exception of P. pinaster which is quite common at sea level (Ruiz de la Torre, 2006). The Third Spanish NFI, conducted between 1997 and 2007, was used to model and predict the site quality in Spanish forests. The NFI plots were only established in woodland areas, according to the FAO definition (FAO, 2001) (i.e. none in non-forested areas) with an intensity of one sampling point every 1 km 2 (1 × 1 km grid). At each sampling point, diameter and height of the trees were measured in four concentric circular subplots with increasing radii from 5 to 25 m. In the 5 m radius subplot, trees with a diameter at breast height (dbh) ≥7.5 cm were measured. In the 10 m radius subplot, trees with dbh ≥12.5 cm were measured. In the 15 m radius subplot, trees with dbh ≥22.5 cm were measured. Finally, in the 25 m radius subplot, trees with dbh ≥42.5 cm were measured. The attributes of the trees measured in each concentric subplot were expanded to per-hectare values by considering a different expansion factor (EF) for each subplot size, which is EF k = 10000/a k where a k is the area of each subplot with size k (k= 5, 10, 15 and 25). Thus, EF k corresponds to the trees per hectare represented by every tree measured in each subplot with size k. (Alberdi Asensio et al., 2010; Hernández et al., 2014). 2.2. Site form models Site quality models were developed for each EFT using the NFI plots with basal area dominance of the main species greater than 90% (Table 1). Site quality was estimated using stand level height-diameter relationships (Huang and Titus, 1993; Stout and Shumway, 1982), which is referred to as site form (Vanclay, 1994) when using dominant diameter (D 0 , cm) and dominant height (H 0 , m). D 0 and H 0 were computed for each NFI plot as the arithmetic mean of the dbh (cm) and tree height (h, m) respectively of the 100 trees per hectare with the mean largest dbh. As we have only one measurement of these variables for each sample plot, the site form models were developed in two steps by using the guide-curve method (Clutter et al., 1983). In the first step a model was fitted to the pairs of D 0 and H 0 , thereby obtaining the average H 0 /D 0 relationship for each EFT. In the second step, the family of site form curves were generated making each of the model parameters dependent on the site form, which is the value of H 0 at a given reference value of D 0 . Thus, each base model generates as many site form family curves as parameters has in the base model formulation. Two well-known base models used in the development of site quality models were considered in the fitting process, namely, the Hossfeld II model (H-II, Eq. (1)) and the Bertalanffy-Richards model (B-R, Eq. (2)), and all possible resultant families of site form family curves were generated for each EFT. The base models were adapted in order to attain predictions of H 0 =1.3 for D 0 = 0. The final model for each EFT was selected according to the performance of the site form curves over the data by considering as criteria the polymorphism of the resulting curves as well as the presence of multiple asymptotes. The reference D 0 was selected to improve the accuracy of predictions, reducing the prediction bias associated with stands where D 0 differs greatly from the reference D 0 (Weiskittel et al., 2011). Furthermore, the range of the diameters observed in the different EFTs and the error of the model in different regions of the validity range of the model were also taken into account. = + + HD a b D 1.3 ( · ) 00 2 02 (1) = +H a b D1.3 ·(1 exp( · ))c 0 0 (2) where H 0 is the dominant height (m), D 0 is the dominant diameter (cm) and a,band care model parameters to be estimated. After selecting the best models for each EFT, we estimated the site form value in the NFI plots. The site form predictions for each EFT were then divided into four site quality classes (actual site quality) using the minimum, the first quartile, the median, the third quartile and the maximum of the site form dataset as breakpoints. Class A represents the highest site quality class, class B and C are the high- and low-inter- mediate site quality classes respectively and, finally, class D is the lowest site quality class. 2.3. Spatial analysis In ecology, as well as in other fields such as forestry, non-linear relationships between variables are common and, therefore, linear modelling often performs poorly (Faraway, 2006). Hence, approaches, such as additive models (Hastie and Tibshirani, 1989; Wood, 2006), which describe complex relationships between the response and the predictors are especially useful. Thus, we fitted the following spatial additive models to predict the form index for each of the five EFT in Spain: = + + + = SFI f X Y f x( , ) ( ) i i i j i J ij i (3) where SFI is the site form index, αrepresents the intercept of the model, f X Y( , ) i i is a spatial smooth function to account for the spatial pattern of the SFI and eliminate the spatial correlation. X i and Y i are the coordinates in meters (UTM, datum ED50 zone 30N) of the i-th plot. f(x ij ) are smooth functions (from j= 1 to J) of physiographic (altitude in m, slope in degrees, and aspect in degrees) and climatic variables to be included in the model. Physiographic variables were extracted from a digital elevation model. We obtained the climatic data from the 200 × 200 m Spanish climatic grid (Gonzalo, 2010). This grid provides both rainfall (total, summer, winter, spring and autumn rainfall) and temperature variables. However, temperature variables are highly Table 4 Model selection and Akaike’s Information Criterion for each European Forests Type. In bold the selected model. Model Model terms P. pinaster F. sylvatica Q. pyrenaica P. nigra P. sylvestris (1) α7766.51 10254.05 14607.96 15852.38 27476.10 (2) (1) + f(X i ,Y i ) 7332.67 9805.29 13951.03 14592.00 25631.34 (3) (2) + f(altitude)7151.53 9720.40 13906.19 14335.19 25333.93 (4) (3) + f(aspect) 7146.83 9699.32 13896.34 14301.22 25312.09 (5) (5) + f(slope) 7146.00 9692.03 13892.92 14302.42 25311.93 D. Moreno-Fernández et al. Forest Ecology and Management 417 (2018) 197–207 200 Fig. 2. Estimated centered smooth functions and 95% confidence intervals of the physiographic variables for the site form spatial additive models of the European Forest Types studied. D. Moreno-Fernández et al. Forest Ecology and Management 417 (2018) 197–207 201 correlated to altitude and were not considered in the analysis. The selection of the variables was carried out according to Akaike’s Information Criterion (AIC) using a forward stepwise procedure setting a reduction of five points of AIC as a significant threshold value. All the variables were entered as smooth functions. Finally, ε i is the error term of the model. We represented the smooth functions using thin plate regression splines (Wood, 2003). These splines keep the basis and the penalty of the full thin plate splines (Duchon, 1977) but the basis is truncated to obtain low rank smoothers. This reduces the computational requirements of the smoothing splines and avoids the problems of the knot placement of the regression splines (Wood, 2003). We plotted semivariograms of the residuals of each model to check if the spatial correlation had been eliminated. Additionally, we created semivariogram envelopes after 99 permutations under the assumption of no spatial correlation (Augustin et al., 2009). Once models were fitted, we classified the form index predicted by the spatial additive model into the four site quality classes (predicted site quality) using the same cutoffs as the site form models. The accuracy percentage was then calculated as the ratio of the plots correctly classified to the total number of plots. 2.4. Site quality maps We used the spatial additive models fitted in the previous section to predict the form index and the site quality in all the plots of the NFI where trees of the target species were present (second column in Table 1). We used the following R packages (R Core Team, 2017) to carry out the statistical analysis: “nls”, “mgcv” (Wood, 2011) and “geoR”. The maps were produced in ArcGis 10.2.2. (ESRI, 2014). 3. Results The parameter estimation for the final site form models and the sitespecific parameters for each EFT are shown in Table 2. H-II was the model selected for all the EFTs with the exception of P. nigra, for which B-R was selected. The site form family of curves represents reasonably well the data for the countrywide conditions covered in the NFI data for each EFT (Fig. 1). The percentage of deviance explained by the spatial additive models ranges from 32.5 to 48.6% (Table 3). The structure of the spatial additive models varies from one EFT to another after the forward stepwise variable selection procedure (Tables 3 and 4). The spatial smooth function (f(X i , Y i )) and the smooth function of altitude (f(altitude)) are included in all the models. In addition, to enhance the performance of the model, the spatial term eliminates the spatial correlation of the residuals (see S1 Fig.). The relationship between the altitude and the site form index varies among the five EFT. In the case of P. pinaster, a species well represented from sea level up to altitudes slightly above 1000 m of elevation, f(altitude) reaches values at lower altitudes and decreases progressively at higher altitudes (Fig. 2), with the largest values of the site form index expected to be observed at lower altitudes. In the case of more montane species (P. sylvestris,P. nigra and F. sylvatica), the 95% confidence intervals of the mean widen at lower elevations, and so the site form shows no significant differences for elevations below than 700–800 m. From that altitude upwards, f(altitude) decreases for these three species, but the confidence interval widens again for high elevations, with a threshold at around 1500 m for beech and 1800 m for P. nigra and P. sylvestris. In the case of Q. pyrenaica, a species well represented at altitudes between 400 and 1600 m, the pattern is different, and f(altitude) shows a maximum plateau for elevations between 500 and 800 m; the confidence intervals also widening for lower and higher elevations. Additionally, we found a significant relationship between f(aspect) and site form index for all the species with the exception of P. pinaster.Fig. 2 indicates that the highest values for the site form index are expected to be found at northern exposures (aspect ranging from 300° to 360° and from 0 to 60°) in F. sylvatica,Q. pyrenaica,P. nigra and P. sylvestris EFT. Even so, the difference between the maximum and minimum values of f(aspect) is greater in the cases of F. sylvatica and P. nigra, whereas aspect is less influential for P. sylvestris and Q. pyrenaica. The variable selection procedure reveals a weak association between site form index and slope in F. sylvatica (p-value = 0.0400). In this regard, f(slope) shows a bell-shaped pattern peaking at 22°. Finally, none of the climatic variables significantly reduced the AIC (less than five points). The largest reductions in AIC for the five EFTs spatial models appear after applying the spatial smoothing. This means that the spatial structure accounts for a larger percentage of variance (Table 4). The second most important variable in AIC reduction is f(altitude). Finally, f (aspect) in F. sylvatica,Q. pyrenaica,P. nigra and P. sylvestris models and f(slope) in F. sylvatica play minor roles. Hence, the estimated centered smooth functions of aspect and slope are quite close to zero and the range of variation of these estimated centered smooth functions are narrower than that of altitude leading to stronger relationships between f(altitude) and site quality (Fig. 2). The accuracy of the site quality ranges from 38.2% in the case of Q. pyrenaica to 47.9% for P. nigra (Table 5). The site form predictions using the spatial additive models mainly classified the plots into intermediate quality site classes (B and C) whereas the number of plots classified into the highest and lowest quality classes (A and D, respectively) was lower. This points to the difficulty of predicting cases of productivity levels well above or below the average. The observed site quality classes are shown in Fig. 3. The predicted site quality maps suggest that most of the stands of the five EFTs studied are located in intermediate site quality areas (Fig. 4 and Table 5). The distribution of the quality classes predicted shows more clearly the gradients that depend mainly on elevations and the aspect. Overall, the correspondence between the predicted and observed maps of site qualities is clear (Figs. 3 and 4), even though the model fails to predict classes A and B, which are in fact observed for areas where other classes are more frequent. Table 5 Confusion matrix for each European Forest Type. Accuracy percentage is shown in brackets. Species Site form model (observed site quality) Spatial additive model (predicted site quality) A B C D Pinus pinaster atlantica (45.4%) A 139 121 44 1 B 37 145 114 10 C 11 109 159 26 D 2 59 133 112 Fagus sylvatica (42.7%) A 135 270 73 1 B 34 274 157 13 C 13 193 220 53 D 2 79 209 189 Quercus pyrenaica (38.2%) A 179 428 113 0 B 50 373 296 3 C 19 255 409 39 D 11 153 417 142 Pinus nigra (47.9%) A 418 309 72 4 B 120 388 253 42 C 31 277 406 89 D 6 117 354 326 Pinus sylvestris (45.2%) A 654 464 161 13 B 197 659 400 36 C 47 468 664 113 D 11 225 696 360 A = high quality site, B = intermediate-high quality site, C = intermediate-low quality site, D = low quality site. D. Moreno-Fernández et al. Forest Ecology and Management 417 (2018) 197–207 202 Fig. 3. Observed site quality of the European Forest Types studied. Darker colors indicate low altitudes whereas light intensities denote high altitudes. D. Moreno-Fernández et al. Forest Ecology and Management 417 (2018) 197–207 203 Fig. 4. Predicted site quality of the European Forest Types studied. Darker colors indicate low altitudes whereas light intensities denote high altitudes. D. Moreno-Fernández et al. Forest Ecology and Management 417 (2018) 197–207 204 4. Discussion We estimated the site quality in the five most important EFTs in Spain expressed as the dominant height and dominant diameter relationship and then created site quality maps using spatial additive models. Other authors found higher percentages of accuracy for site index in Mediterranean species: from 64% to 71% in P. sylvestris (Bravo and Montero, 2001; Bueis et al., 2016) and from 61 to 75% in Pinus pinea L. (Bravo-Oviedo and Montero, 2005; Bravo et al., 2011). However, these studies were performed at finer scales, local or regional, whereas our aim was to estimate the site quality at a national scale. In this regard, it has been shown that the relationships between environmental variables and site index are stronger at finer scales (Chen et al., 2002). Our results are in concordance with the autoecology of the species, confirming the suitability of our approach from a biological point of view. The negative relationship between P. pinaster site quality and altitude has previously been reported in Mediterraean areas (Bravo- Oviedo et al., 2011) in other studies under Atlantic conditions (Álvarez- Álvarez et al., 2011; Eimil-Fraga et al., 2014). This species is distributed across areas of Atlantic climate with its central habitat in terms of elevation ranging from 0 to 1000 m, reaching the maximum levels of site productivity at 513 m asl (Álvarez-Álvarez et al., 2011). Therefore is clearly differentiated from the other two pine species studied, which are genuinely montane (Gandullo and Sánchez-Palomares, 1994). In distribution models, altitude or other related variables such as temperature usually show a positive or inverse U-shaped relationship with the occurrence of montane species (Hernández et al., 2014; Moreno- Fernández et al., 2016). That was observed for site productivity in P. rubens (Chen et al., 2002) and, in our study case, in Q. pyrenaica. This species is known to grow at a range of elevations between 400 and 1400 m, with its optimum habitat in northern areas at lower elevations than in the southern ones (Díaz-Maroto et al., 2007; Sánchez-Palomares et al., 2008). An inverse U-shaped relationship may have been expected for the other three species, the central habitats of which are also within elevation ranges classifiable as montane (Gandullo and Sánchez- Palomares, 1994). The occurrence of F. sylvatica at low elevations in Spain is mainly restricted to humid sites on the slopes of northern Mountain ranges, where the species is not restricted by droughts and shows high productivity rates at an elevation range of 400–800 m (Sánchez-Palomares et al., 2004). Some stands showing poor productivity at elevations below the central habitat of these species and located in dryer areas present a mixed composition and therefore, as we only used plots in which the main species of the EFT was dominant, these stands were not considered in the study. Thus, in the plots used in this study, the species grow under conditions which are sufficiently favorable to form monospecific stands or have been favored by forest management. Certain factors associated with systems at high elevations, such as shorter growing periods (Beniston, 2003), may explain the general negative association. The results obtained as regards the influence of aspect on the site form provide new insights into the dependence of stand productivity on site parameters, in addition to those reported in studies of species distribution or regeneration niches (Gómez-Aparicio et al., 2006; Hernández et al., 2014). A clear effect of aspect on stand productivity was particularly evident in the case of F. sylvatica, a shade tolerant species that grows better in northern exposures with more humid conditions (Ruiz de la Torre, 2006). Hernández et al. (2014) reported an increase in the presence of F. sylvatica on north-facing exposures in northern Spain over the last 40 years. Such a clear positive effect of northern exposures on site form was also observed for P. nigra, a species that is considered to have an intermediate shade tolerance, at least in the case of the subspecies salzmannii (Ruiz de la Torre, 2006), although regeneration in this species has been found to be more tolerant of high levels of irradiances than several broadleaved species (Gómez-Aparicio et al., 2006). The role of aspect on site quality was found to be smaller for Q. pyrenaica and P. sylvestris, even though the site form at northern exposures was frequently found to be significantly greater. Cañellas et al. (2000) have already highlighted the fact that the influence of aspect on site index in P. sylvestris site index is expected to be greater for stands with Mediterranean conditions. P. pinaster is a light-demanding species (Ruiz de la Torre, 2006) and the subspecies atlantica requires humid conditions for a good development. It grows, as mentioned above, mainly at relatively low altitudes close to the coast (Álvarez- González et al., 2005), influenced by the Atlantic climate with mild summer temperatures (Elena Roselló, 1997); therefore, it is not surprising that aspect had no statistically significant influence on site quality in this species. Other local variables, such as soil variables, are expected to influence the site index (Álvarez-Álvarez et al., 2011; Bravo et al., 2011; Bueis et al., 2016; Chen et al., 2002). Soil variables, however, show great spatial variations (see Vanderlinden et al. (2003) for soil water holding in southern Spain). Gathering soil data, such as soil texture, soil depth, organic matter or nutrients content, for the whole country would be very useful for predicting site quality but it would be highly timeconsuming and expensive (McBratney et al., 2006). The unmeasured variables could in fact explain the lower ability of the model to classify correctly some of the A and D productivity class plots. Nevertheless, the effect of unmeasured environmental variables was partially accounted for by the spatial smooth function (Moreno-Fernández et al., 2018). Additionally, the spatial smooth function eliminated the spatial correlation minimizing the type I error rates (Dormann et al., 2007). Furthermore, if we consider the strong rainfall gradients in Spain (Muñoz- Díaz and Rodrigo, 2004), the lack of significant relationships between rainfall variables and site quality would be unexpected. However, physiographical variables and the spatial smooth function may, in fact, be absorbing the effect of the rainfall. It should also be mentioned that both F. sylvatica and Q. pyrenaica are able to sprout and thus can form coppice forest (Ruiz de la Torre, 2006). Hence, coppicing management systems have traditionally been those most frequently employed in Spanish stands of Q. pyrenaica (Adame et al., 2008). It is known that the growth pattern of a given species differs from coppice to high forest, with the former showing faster initial height growth although this growth is less sustained than in high forest (Ciancio et al., 2006; Haneca et al., 2005). However, the effect of stand structure on site form is not clear. Further research is needed to address this aspect and if deemed necessary, stand structure can be integrated into site quality models (Adame et al., 2008). Furthermore, within the distribution area of each species in Spain, several different subspecies can be found. This is particularly evident in the case of P. nigra which presents high intraspecific variability within natural stands (Ruiz de la Torre, 2006). Moreover, the use of Corsican and Austrian provenances in reforestation has resulted in still further variability. This can lead to the attribution of different site qualities for each subspecies at the same location when considering dominant height and age as site quality criteria (Moreno-Fernández et al., 2014). However, since the fact that the P. nigra subspecies appear rarely in close proximity to each other, spatial smoothing can capture the influence of the subspecies on site form. In addition, the use of the dominant diameter instead of age can absorb part of the effect of the subspecies differentiation if the subspecies only differ in growth rates. Our results provide a tool which could be used by forest managers in land use planning as well as in forest policy decision-making at a national scale. For instance, the proposed site-form models may be used for estimating site quality in even- or uneven-aged stands of the studied species, and the maps for large-scale forest planning. Additionally, this methodology fulfills the current demands for reporting forest information according to international criteria and provides a viable approach for standardizing the method used to estimate site quality, thereby producing internationally comparable data. This methodology could be used in other countries and the maps could be extended to the whole of Europe. Thus, it would be possible to assess the way in which D. Moreno-Fernández et al. Forest Ecology and Management 417 (2018) 197–207 205