DEVELOPING AN INDEX FOR FOREST PRODUCTIVITY MAPPING - A CASE STUDY FOR MARITIME PINE PRODUCTION REGULATION IN PORTUGAL

DESENVOLVIMENTO DE UM ÍNDICE PARA O MAPEAMENTO DA PRODUTIVIDADE FLORESTAL - UM ESTUDO DE CASO PARA A REGULAÇÃO DA PRODUÇÃO DO PINHEIRO BRAVO

ABSTRACT

Productivity is very dependent on the environmental and biotic factors present at the site where the forest species of interest is present. Forest site productivity is usually assessed using empirical models applied to inventory data providing discrete predictions. While the use of GIS-based models enables building a site productivity distribution map. Therefore, the aim of this study was to derive a productivity index using multivariate statistics and coupled GIS-geostatistics to obtain a forest productivity map. To that end, a study area vastly covered by naturally regenerated forests of maritime pine in central Portugal was used. First, a productivity index (PI) was built based on Factorial Correspondence Analysis (FCA) by incorporating a classical site index for the species and region (Sh25 - height index model) and GIS-derived environmental variables (slope and aspect). After, the PI map was obtained by multi-Gaussian kriging and used as a GIS layer to evaluate maritime pine areas by productivity class (e.g., low, intermediate and high). In the end, the area control method was applied to assess the size and the number of compartments to establish by productivity class. The management compartments of equal productivity were digitized as GIS layer and organized in a temporal progression of stands’ age regularly available for cutting each year during a 50-year schedule. The methodological approach developed in this study proved that can be used to build forest productivity maps which are crucial tools to support forest production regulation.

Keywords:
Site Productivity Classes; Area Control Method; Forest Management Compartments

RESUMO

A produtividade florestal depende dos fatores ambientais e bióticos da estação onde a espécie florestal se encontra. A produtividade da estação florestal é geralmente avaliada com modelos empíricos aplicados a dados de inventário fornecendo previsões discretas. Enquanto o uso de modelos num SIG permite a construção de mapas de distribuição de produtividade da estação. Assim, o objetivo deste estudo foi desenvolver um índice de produtividade utilizando estatística multivariada e geoestatística acopladas ao SIG para obter um mapa de produtividade florestal. Para o efeito, selecionou-se uma área de estudo com floresta de regeneração natural de pinheiro bravo no centro de Portugal. Primeiramente, o índice de produtividade (PI) foi construído usando a Análise Fatorial das Correspondências (FCA) com um índice de qualidade de estação para a espécie e região (Sh25- modelo de índice de altura) e variáveis ambientais (declive e exposição). Depois, o mapa PI foi obtido por Krigagem multi-Gaussiana e usado como camada SIG para avaliar as áreas da espécie por classe de produtividade (p.e., baixa, intermedia e alta). No final, o método de controle da área foi aplicado para determinar o tamanho e o número de parcelas a serem estabelecidas por classe de produtividade. As parcelas de igual produtividade foram digitalizadas em SIG e organizadas segundo uma progressão de idades disponíveis regularmente para abate anual ao longo de 50 anos. A abordagem metodológica desenvolvida neste estudo provou ser útil na construção de mapas de produtividade florestal que são instrumentos fundamentais para apoiar a regulação da produção florestal.

Palavras-Chave:
Classes de Produtividade da Estação; Método do Controlo da Área; Parcelas de Gestão Florestal

1.INTRODUCTION

Knowing the site productivity of a forest species in a specific area of interest is essential to attain sustainable forest management (Mitsuda, 2007; Cruz et al., 2008; Ribeiro et al., 2016; Yue et al., 2016; Dãnescu et al., 2017). Therefore, the practical need for forest productivity assessment has encouraged the development of several quantitative methodologies. Forest site productivity is usually measured indirectly as a function of site characteristics and/or vegetation as it depends both on natural factors inherent to the site (e.g. topographic, climatic and soil factors) and on management-related factors. The most widely used and accepted site productivity indicator is the site index (Vanclay, 1994; Skovsgaard and Vanclay, 2008; Bontemps and Bouriaud, 2014). However, site index does not represent the true site quality as it does not fully integrate the site factors. Therefore, more recently prediction models for site index evaluation using regression techniques have included environmental factors such as climatic, topographic or soil related attributes to improve forest productivity assessment (Fontes et al., 2004; Aersten et al., 2010; Salas, 2010; Nunes et al., 2011; Dãnescu et al., 2017).

Relating site productivity to environmental factors, derived from digital terrain analysis using GIS-based models, can easily provide site productivity distribution maps for wide areas, which are useful and powerful tools for sustainable forest management as well as for landscape management (Mitsuda et al., 2007). Current prediction models are applied to inventory data and, as a result, only provide a punctual prediction of site productivity. The application of geostatistical techniques to point-data GIS layers enables, in one hand, to move from a punctual representation of forest attributes to a continuous representation as interpolated surfaces and, in another hand, to obtain a measure of uncertainty attached to the prediction of unsampled sites (Vanclay, 1994; Goovaerts, 1997; Santos et al., 2003; Pelissari et al., 2015). Coupled GIS and geostatistical methodologies have been successfully used to predict several forest variables, namely for forest productivity mapping (Hock et al., 1993; Payn et al., 1999; Bognola et al., 2009). At national or regional level studies, accurate kriging results are easily obtained as a wide range of biophysical conditions (e.g. soil, climate and topographic features) are in consideration (Mandallaz, 2000; Nanos et al., 2003; Lima et al., 2006; Benavides et al., 2009; Bognola et al., 2009; Hlásny et al., 2017). However, for local studies, a microclimatic approach should be explored instead to incorporate variables such as aspect and factors contributing greatly to the development of the soil, such as slope and elevation, especially in mountainous regions (Benavides et al., 2009).

The species maritime pine (Pinus pinaster Aiton) has a great expression in Portuguese forest area (714 10³ ha; 23%) (ICNF, 2013) and is mainly located at the center and north of the country. It is a Mediterranean species that thrives well until 700-800 m preferably in Atlantic influence facing slopes(Alveset al., 2012). In Portugal, several classical site index models are available to support this species management (Páscoa, 1990; Nunes et al., 2011; Alegria and Tomé, 2013). Regarding Portuguese forest management plans mapping management compartments for forest production regulation is necessary. The area control method (Davis and Johnson, 1987) is the easiest way to regulate an unmanaged forest and guarantees that the regulated structure is attained within one rotation (i.e. the number of years to final harvest). This method has also the virtue of portraying the harvest in terms of area to cut. Since forest land always displays some site productivity variability, areas of equal productivity rather than of equal surface must be evaluated.

Therefore, this paper presents a new approach to derive a productivity index using multivariate statistics and coupled GIS-geostatistics, as opposed to previous studies that used regression techniques, to obtain a forest productivity map to support forest production regulation. In that sense, the first hypothesis addressed in the study was that a productivity index (PI) could be derived by incorporating a classical site index and topographic factors (slope and aspect) through the application of Factorial Correspondence Analysis (FCA). The second hypothesis was that a prediction map of the PI could be built using coupled GIS-geostatistics analysis. Finally, the third hypothesis was that the PI map could be used as a GIS layer for forest productivity assessment to support forest regulation. To test these hypotheses, a study area vastly covered by naturally regenerated forests of maritime pine in the center region of Portugal was selected as a case study.

2. MATERIALS AND METHODS

2.1. Study area and data

The study area is a 3,100 ha municipality located in central Portugal (Fig. 1a). The mountain area at the West (Serra do Moradal) is included at the Geopark of “Naturtejo da Meseta Meridional” which is integrated in UNESCO World Geoparks since 2006. Most of this mountain area is covered by forests and semi-natural areas (90%), where forests are mainly composed of pure maritime pine (78%) (DGT, 2017). The northern zone of the study area was burnt in 2003 (ICNF, 2017) and is now covered by young natural regeneration of maritime pine open forest or scrubland whereas the southern zone is mainly occupied by mature maritime pine forest (Fig. 1b). A geo-referenced 500 m spacing grid was laid over to the study area in order to locate point-support samples (one sample plot every 25 ha) in the maritime pine forest (Fig. 1d). Due to the fire in 2003, no inventory data was available at the 31 plots in the northern zone. These plots were either of young naturally regenerated maritime pine stands that were too small to be measured (18 plots) or of scrubland (13 plots). In the southern zone, covered with mature maritime pine, 50 plots each of 500 m² were installed for inventory data collection in 2010. A total of 2,259 trees were sampled in these plots and their diameters at breast height (d), i.e. 1.30 m above ground, were measured. Total heights (h) were measured for a subset of 497 trees selected proportionally to the class diameter frequency in each sample plot. The data collected allowed the evaluation of several variables that describe stands growth conditions, site productivity and wood production.

Figure 1
Study area and input data: a) geographic location and maritime pine distribution in Portugal (Godinho-Ferreira et al., 2005); b) maritime pine forest and open forest in 2007 and burnt area in 2003; c) digital elevation model (DEM) - spatial resolution 1 m; d) inventory data - Sh25 classes; and GIS-derived topographic data - slope classes and aspect classes.
Figura 1
Área de estudo e dados de entrada: a) localização geográfica e distribuição do pinheiro bravo em Portugal (Godinho-Ferreira et al., 2005); b) florestas de pinheiro bravo em 2007 e área ardida em 2003; c) modelo digital de terreno (MDT) - resolução espacial de 1 m; d) dados de inventário - classes de Sh25; e variáveis topográficas obtidas em SIG - classes de declive e classes de exposição das encostas.

These maritime pine stands are uneven-aged and so site productivity was assessed using a height index model (Sh25) developed for uneven-aged stands of the species and region(Alegria and Tomé, 2013). This height index model (Clutter et al., 1983; Huang and Titus, 1993; Vanclay, 1994) for productivity evaluation was adjusted by regression techniques using the following monomolecular function modified by Meyer (1940):

(Eq.1) $\mathit{h}=1.3+\mathit{A}\left(1-{\mathit{e}}^{-\mathit{kd}}\right)$

with, A - asymptote and k- species coefficient. The height index curves were then obtained using the guide curve method as described in Clutter et al. (1983):

(Eq2) $\mathit{Sh}25=1.3+\left(\mathit{h}-1.3\right)\frac{\left(1-{\mathit{e}}^{-1.1725}\right)}{\left(1-{\mathit{e}}^{-0.0469\mathit{d}}\right)}$

where, h is the tree total height (m) and d is the tree diameter, over bark, at breast height (cm). The Sh25 index, also called site form by Vanclay (1994), is defined as the expected tree height at the reference diameter of 25 cm and was computed as the average of individual tree values. This Sh25 index is positively correlated with both the mean annual increment of stand volume (MAI) and the classical site index for evenaged stands of the species (SI50) (Páscoa, 1990). However, the Sh25 index has the advantage of not requiring age data by contrast to the site index SI50.

The topographic variables considered for this local approach were the slope and aspect as suggested by Benavides et al. (2009). A digital elevation model (DEM) with a resolution of 1 m was used to compute both slope and aspect layers using the 3D Analyst tool in software ArcGIS Desktop v9.3 (ESRI, 2004) (Fig. 1c).

2.2 Productivity index - multivariate statistics approach

The Sh25 index (i.e. a plant-based variable, evaluated at a fine scale but without continuous spatial coverage) and the topographic variables slope and aspect (i.e. GIS digital data of continuous spatial coverage) were used to develop a productivity index (PI). The topographic variables slope and aspect were GIS-derived for the sample plots. The three variables were first reclassified in three classes: (1) low, (2) intermediate and (3) high. Sample plots’ Sh25 index was coded as: (1) Sh25<15 m, (2) Sh25Î [15, 18[m and (3) Sh25e”18m (Fig. 1d). Slope is indirectly related to soil features such as vertical characteristics, external drainage, water storage, effective thickness, and rocky outcrops. Therefore, sample plots’ slope was coded as: (1) slope Î [10-15[%, (2) slope Î [15-35[% and (3) slope e” 35% (Fig. 1d). It is commonly accepted that northern facing slopes are colder and southern facing slopes are hotter which has an impact on species distribution and development. So, sample plots’ aspect was coded as: (1) slopes facing NW-N-NE (orientation between 0º- 67.5º and 292.5º-360º), (2) slopes facing E-SE (orientation between 67.5º-157.5º) and (3) slopes facing S-SW-W (orientation between 157.5º-292.5º) (Fig. 1d). Slope class and aspect class for each sample plot (500 m² of the area) was obtained by a majority filter applied for the cells occupying each plot (zonal stats function in ArcGIS (GRID).

The productivity index (PI) was built using a two- step process based on Factorial Correspondence Analysis (FCA). In a first step, data was assembled in a complete disjunctive matrix, which allocates for each attribute each sample plot into three classes corresponding to the lower (1), intermediate (2) and higher values (3). The main advantage of FCA is that symmetry is conferred to the data matrix, thus allowing the simultaneous study of correlations within and between variables and samples and a hierarchical sample interpretation according to their level is possible (Pereira et al., 1993, 2015). This procedure showed that low Sh25 (Sh25<15 m) is associated with high slope (Slope e”35%) and slopes facing E-SE (orientation between 67.5º-157.5º) whereas high Sh25 (Sh25e”18 m) is associated with low slope (Slope Î [10-15[%) and slopes facing S-SW-W(orientation between 157.5º-292.5º) (Fig. 1d).

In a second step, based on the previous attributes relationships analysis, the productivity index (PI) was built using FCA tailored as a discriminant procedure. The method proceeds as follows: (1) selection of the parameters to be included in the index; (2) standardization of the parameters and (3) aggregation of the parameters. The aggregation of the standardized parameters was developed by Benzécri in the early seventies (Benzécri, 1977): The Factorial Correspondence Analysis (FCA) belongs to a group of factor extraction methods whose main objective is to discover the underlying pattern of relationships within a data set. This is basically done by rearranging the data into a smaller number of uncorrelated ‘‘components” or ‘‘factors” that are extracted from the data by statistical transformations. Such transformations involve the diagonalization of a similarity matrix of the variables, such as a correlation or variance-covariance matrix. Each factor describes a certain amount of the statistical variance of the analyzed data and is interpreted according to the inter-correlated variables. A detailed discussion of the theory behind FCA goes beyond the scope of this article, but its application in the present case is rather straightforward.

The method used here (Pereira et al., 1993; Albuquerque et al., 2010; Pereira et al., 2015) subdivides each of the three target variables (Sh25, slope and aspect) into three classes using the procedure. Each sample is located on an arbitrary dimensionless scale defined by two extreme poles (virtual vectors) one of very high quality and the other one of very low quality. Following this procedure, each experimental sample result is located on an arbitrary scale deûned by these two extreme poles using the FCA supplementary projection. The resulting scores correspond to the ûnal index values which range between ‘-1’ to ‘1’ (low to high). The new variable (PI) obtained by this method is defined as a regionalized variable and it is additive by construction: the barycentric principle of FCA guarantees that two samples with given profiles in the variables can be replaced by a new individual, whose profile is given by the coordinates of the center of gravity of the two original samples (Benzécri, 1982). The developed productivity index PI ranks the combination of the variables Sh25, slope and aspect.

2.3 Mapping productivity index - geostatistics approach

The productivity index (PI) was mapped using multi-Gaussian kriging implemented in SpaceStat v2.2.17 (BioMedware, 2011) and that proceeds in three steps: 1) normal score transforms of the PI data, 2) interpolation of normal scores using ordinary kriging, and 3) back- transform of the results using the empirical procedure developed by Saito and Goovaerts (2000). The PI values were transformed into normal scores to attenuate the impact of extreme values on the computation of the variogram). Under the multi-Gaussian model, each Conditional Cumulative Distribution Function (CCDF) is Gaussian and has for mean and variance, respectively the kriging estimate and kriging variance. The quality of the model was assessed using a cross-validation or leave-one-out approach that proceeds as follows: 1) each PI observation is removed at a time and the CCDF is modeled using neighboring data, 2) p-probability intervals, centered around the median of the CCDF are then built; for example, 0.5 is bounded by the lower and upper quartiles of the CCDF. The scattergram of the estimated versus expected proportions is called, “accuracy plots” which were here built using the program ACCPLT from GSLIB (Deutsch and Journel, 1992). The closeness of the observed and expected proportions was quantified using the goodness statistic (best value =1) (Goovaerts, 2001) that penalizes more the situation where the observed proportion is lower than expected (inaccurate case).

2.4 Forest production regulation - management compartments

To achieve a fully-regulated forest for maritime pine in the southern zone the area control method was applied. A 50-year rotation was considered taking in consideration the sylvicultural prescription by Oliveira (1999), which was proved to be the most appropriate for this forest (Alegria, 2011), to attain a sustainable forest management. This sylvicultural prescription proposes the final harvest at 45 years. However, five more years were considered to guarantee the species natural regeneration after harvest resulting in the 50- year rotation abovementioned. Using this sylvicultural prescription it is expected that a mean annual increment of stand volume at the rotation age of 45 years, of 5 m³ha^-1yr^-1 in low quality sites (Sh25<15), 8 m³ha^-1yr^-1 in intermediate quality sites (Sh25Î [15, 18[ m) and 12 m³ha^-1yr^-1 in high-quality sites (Sh25e”18) respectively, will be observed (Alegria, 2011).

The produced PI map was used as an input GIS layer to obtain maritime pine areas by productivity class, in the southern zone, the PI map was reclassified into three classes: low (1) for PI<-0.48, intermediate (2) for PIÎ [-0.48, -0.16[and high (3) for PIe”-0.16. This classification was possible due to the high association between the PI and the Sh25. Then, the maritime pine areas by PI class were evaluated and subsequently used under the area control method to assess the number of hectares to harvest per year (i.e. the compartments’ size) and the years of harvesting (i.e. the compartments’ number) for each productivity class (Table 1). This resulted in a total of 50 compartments: four compartments of 43 ha in low PI class, 18 compartments of 28 ha in intermediate PI class and 28 compartments of 19 ha in high PI class. This information supported the vectorization of the layer management compartments.

Table 1
Area control method (Davis and Johnson, 1987) for maritime pine using PI for site productivity assessment and a 50-year rotation sylvicultural prescription.
Tabela 1
Método do controlo da área (Davis And Johnson,1987) para o pinheiro bravo utilizando o PI para determinar a produtividade e uma prescrição silvícola de revolução de 50 anos.

Site	Area of forest	MAI45(m3ha-1yr-1)	Relative productivity	Area of equal productivity	Long term sustainability	Productivity equivalent*	Hectares to harvest per year	Years of harvesting
PI class	(ha)			(ha)	(m3)	(ha)
(1)	(2)	(3)	(4)=(3)/8	(5)=1/(4)	(6)=(2)x(3)	(7)=(2)x(4)	(8)= s x(5)	(9)=(2)/(8)
Low	172.6	5	0.6	1.6	863	107.9	45.0	3.8
Intermediate	510.4	8	1.0	1.0	4083	510.4	28.1*	18.1
High	525.8	12	1.5	0.7	6310	788.8	18.8	28.0
Total	1,208.8				1,1256	1,407.0	*s= 1,407/50	50

• Legend: (1) Site productivity PI class – is classified as low for PI<-0.48; intermediate for PI Î [-0.48,-0.16[and high for PIe"-0.16; (2) Area of forest (ha) – is the maritime pine area by PI class evaluated using the PI map; (3) ₄₅ – is the mean annual increment of stand volume at the rotation age of 45 years (m³ha^-1yr^-1) assigned to each PI class; (4) Relative productivity – is obtained as the ratio between "₄₅" in each PI class and the intermediate value of 8 m³ha^-1yr^-1; (5) Area of equal productivity – is obtained as the inverse of "Relative productivity" values; (6) Long term sustainability – is obtained as the product of the "Area of forest" with the "₄₅" in each PI class; (7) Productivity equivalent – is obtained as the product of the "Area of forest" with the "Relative productivity" in each PI class; (8) Hectares to harvest per year – is obtained as the product of the "Area of equal productivity" with the equal productive area (i.e.,s= area of forest divided by the rotation age); (9) Years of harvesting – is obtained as the ratio between the "Area of forest" and the "Hectares to harvest per year" in each PI class.

The management compartments’ boundaries were drawn over the maritime pine areas in the southern zone taking in consideration the network of roads, rivers, and streams. Compartments greater than 20 ha were divided into two or more sub-compartments to mitigate the risk of erosion after harvest. Last, the compartments were organized in a temporal management schedule of 50 years to achieve a fully-regulated forest within this rotation.

Stand’s age in each compartment was obtained using the nearest sample plot data as a reference, knowing however that sacrifices are always needed to be made, that there will be stands either younger or older than 45 years that will be harvested until a fully-regulated forest is attained.

Finally, to set the harvest schedule for the next 50 years the following harvest rules were used: 1. the rotation age is 50 years; 2. the area regulation is attained by cutting and regenerating one compartment of equal productivity each year over the 50 year rotation; 3. the compartment to harvest is always the one with the highest age; 4. harvesting adjacent compartments is avoided and a dispersed structure of harvested compartments is promoted to mitigate erosion risk and other environmental impacts; 5. compartments are cut and then grow over the 50 year rotation.

3. RESULTS

Factorial Correspondence Analysis (FCA) allowed developing the productivity index PI by incorporating the Sh25 index, slope and aspect for maritime pine productivity assessment. Since the data were collected on a 500 m spacing grid, the experimental variogram of PI data could not be computed for distances smaller than 500 m. However, the availability of individual tree statistics for each plot allowed the estimation of the within-site variability, which was used as an estimate of the nugget effect. The final model thus consisted of a spherical model with a range of 1 km and a small nugget effect (less than 10%). This model was used to interpolate PI values to the nodes of a 500 m spacing grid covering the southern zone of the study area.

The resulting map was produced with a resolution of 70 m (Fig. 2a) and clearly shows low productivities in the south-western zone and high productivities in the northern zone of the study area. The quality of the model of uncertainty created using multi-Gaussian kriging was assessed using accuracy plots created by a cross-validation approach. The results of cross- validation statistics obtained for the forest quality index (PI) showed a good agreement between observed and expected proportions of observations within PIs, which was confirmed by the high value of the goodness statistic (0.857).

Figure 2
Maritime pine in the southern zone of the study area: a) productivity index (PI) map created by multi-Gaussian kriging; b) PI classes map; c) management compartments map by PI class (i.e.,50 compartments of equal productivity); d) harvesting schedule to attain a fully-regulated forest within a 50-year rotation (by decades).
Figura 2
Pinheiro bravo na zona sul da área de estudo: a) mapa da distribuição espacial do índice de produtividade obtido por krigagem multi-Gaussiana; b) mapa das classes de PI; c) mapa das parcelas de gestão por classes de PI (i.e.,50 parcelas de igual produtividade; d) Plano de cortes para a regulação da produção florestal dentro de uma revolução de 50 anos (por década).

The PI map (Fig. 2a) after being reclassified in three productivity class (low - PI<-0.48; intermediate - PI Î [-0.48, -0.16[; and high - PIe”-0.16) enabled the evaluation of maritime pine areas by productivity class (Fig. 2b). Maritime pine area in the southern zone (1,209 ha) was divided in 50 compartments of equal productivity: four compartments of 43 ha located in low PI class (173 ha), 18 compartments of 28 ha located in intermediate PI class (510 ha) and 28 compartments of 19 ha located in high PI class (526 ha) (Fig. 2c). The management compartments have an irregular shape because their boundaries were drawn taking in consideration the network of roads, rivers, and streams. Moreover, most of the compartments were split into patches smaller than 20 ha and set apart from one another to reduce the risk of post-harvest erosion (Fig. 2d). In the end, the compartments were organized in a temporal management schedule of 50 years by decades (Fig. 2d) in respect to the harvest rules mentioned earlier, so that a fully-regulated forest will be attained within one rotation of 50 years. The first compartment scheduled to be harvest is compartment number 16. In the following year, compartment number 37 will be harvest, and this process continued for all compartments during the 50 year rotation. As a result, in the first decade the sequence of compartments to be harvest are as follows: 16, 37, 4, 5, 6, 39, 20, 47, 19 and 48. In the second decade are the following: 49, 42, 40, 24, 44, 28, 25, 45, 36 and 35. In the third decade, the following: 38, 31, 2, 3, 27, 9, 23, 21, 22 and 41. In the fourth decade, the following: 33, 30, 7, 26, 34, 46, 11, 13, 1 and 18. And in the fifth decade, the following: 17, 8, 32, 12, 10, 29, 15, 43, 50 and 14 (Fig. 2c; Fig. 2d). After this 50-year rotation, this forest is fully-regulated and a new cycle can be re-started.

4. DISCUSSION

A productivity map for maritime pine was created using a productivity index (PI) to support forest regulation and to attain sustainable management. This index incorporates the Sh25 index and the topographic variables slope and aspect. Multi-Gaussian kriging (interpolation) was used for the PI map construction, over the southern zone of the study area. The estimated map allows productivity prediction for the species which are consistent and biologically realistic. The species thrives very well in all study area but may have limitations above an elevation of 800 m. and prefers Atlantic influence facing slopes (Oliveira et al., 2000; Alves et al., 2012).

That preference was identified by the Factorial Correspondence Analysis (FCA) in which the highest values of Sh25 index (Sh25e”18 m) were associated with slopes facing S-SW-W (orientation between 157.5º- 292.5º). Moreover, climate (temperature and rainfall) and soil limiting factors for maritime pine development are very similar over the study area. Therefore, the GIS-derived topographic variables slope and aspect along with stand condition seem to be the main culprits for the observed differences in productivity. These results are in accordance with the use of a microclimatic approach for local studies as suggested by Benavides et al. (2009).

The PI map used as GIS layer proved to be suitable for landscape planning and forest management support. Both compartments’ size and number by PI class were obtained through the application of the area control method (Davis and Johnson, 1987). This allowed organizing these compartments of equal productivity in a temporal management schedule of 50 years so a fully-regulated forest will be attained after one rotation. Additionally, such a progression of stands’ age classes changes the spatial and temporal pattern of the forest which will offer an avenue to improve landscape diversity (Kerr, 1999; Grant et al., 2012). Compartments have an irregular shape, are composed of sub-compartments smaller than 20 ha and set apart from one another to reduce erosion risk after harvest and to improve landscape diversity. In fact, the application of the method of the area control offers the opportunity to consider species conversion, after each compartment being harvest, for landscape fragmentation by mixing compartments of different combustibility. Therefore, maritime pine forest can be converted to a fragmented landscape structure by introducing the recommended native broadleaves species for this region. As a result, other than fire hazard mitigation related benefits could be provided, namely landscape diversity and scenery (Navalho et al., 2017).

While implementing this regulation method, it is advisable to monitor the timber production (e.g. each five-ten-year period) so the proposed site productivity classes (PI classes - MAI₄₅) can be validated and more detailed and accurate PI classes may be defined in the future.

The PI map obtained in this study is an important support tool in Portuguese Forest Management Plans to derive the final map of management compartments for forest regulation and to attain sustainable forest management. This management compartments’ layer will respond to the management questions of “how much”, “where” and “when” to harvest. The sylvicultural prescription by Oliveira (1999) will respond to the management question of “how” to harvest. Finally, since study area’s maritime pine forest is mainly private, to respond to the last management question of “who” is going to harvest, creating a Forest Intervention Zone gathering the forest owners into this common management plan would be advisable.

5. CONCLUSION

The methodology developed in this study stands as an alternative approach to the classical site index modeling and has the advantage of using an available classical index for the species which is then improved by incorporating GIS-derived topographic variables (slope and aspect) using Factorial Correspondence Analysis. This methodological approach can be used to produce forest productivity maps by coupled GIS- geostatistics to assist decision-making at a regional level on landscape planning for sustainable forest management.

6. REFERENCES

• Aersten W, Kint V, van Orshoven J, Özkan K, Muys, B. Comparison and ranking of different modeling techniques for prediction of site index in Mediterranean mountain forests. Ecological Modeling. 2010;221(8):1119-30.
• Albuquerque T, Dias VH, Poellinger N, Pinto JP. Construction of a quality index for granules produced by fluidized bed technology and application of the correspondence analysis as a discriminant procedure. European Journal of Pharmaceutics and Biopharmaceutics. 2010;75(3):418-24.
• Alegria C. Modelling merchantable volumes for uneven aged maritime pine (Pinus pinaster Aiton) stands established by natural regeneration in the central Portugal. Annals of Forest Research. 2011;54(2):197-214.
• Alegria C, Tomé M. A tree distance-dependent growth and yield model for naturally regenerated pure uneven-aged maritime pine stands in the central inland of Portugal. Annals of Forest Science. 2013;70:261-76.
• Alves AM, Pereira JS, Correia AV. Silvicultura. A gestão dos ecossistemas florestais. Lisboa: Lisboa: Fundação Calouste Gulbenkian; 2012. 597p.
• Benavides R, Roig S, Osoro K. Potential productivity of forested areas based on a biophysical model. A case study of a mountainous region in northern Spain. Annals of Forest Science. 2009;66:108.
• Benzécri JP. L’analyse des correspondances. Les Cahiers de l’Analyse des Données. 1997;2(2):125-42.
• Biomedware. Space Stat User Manual version 2.2. BioMedware; 2011. 320p
• Bognola IA, Ribeiro PJ, Silva EAA, Lingnau C, Higa AR. Modelação uni e bivariada da variabilidade espacial do rendimento de Pinus taeda L. Floresta. 2008;38(2):373-85.
• Bontemps JD, Bouriaud O. Predictive approaches to forest site productivity: Recent trends,challenges and future perspectives. Forestry. 2014;87(1):109-28.
• Clutter JL, Fortson JC, Pienaar LV, Brister GH, Bailey RL. Timber management: a quantitative approach. New York: John Wiley & Sons; 1983. 307p.
• Cruz JP, Leite HG, Soares CPB, Campos JCC, Smit L, Nogueira GS. Curvas de crescimento e de índice de local para povoamentos de Tectona grandis em Tangará da Serra, Mato Grosso. Revista Árvore. 2008;32(4):679-85.
• Dãnescu A, Albrecht AT, Bauhus J, Kohnle U. Geocentric alternatives to site index for modeling tree increment in uneven-aged mixed stands. Forest Ecology and Management. 2017;392:1-12.
• Davis L, Johnson K. Forest management. New York: McGraw-Hill, 1987. 256p.
• Deutsch CV, Journel AG. GSLIB: Geostatistical Software Library and User’s Guide. London: Oxford University Press; 1992. 369p
• DGT. Carta de ocupação do solo 2007 para Portugal. [acessado em:05 jun. 2017] Disponível em: http://www.dgterritorio.pt/cartografia_e_geodesia/cartografia/cartografia_tematica/carta_de_ocupacao_do_solo__cos_/cos__2007/
  » http://www.dgterritorio.pt/cartografia_e_geodesia/cartografia/cartografia_tematica/carta_de_ocupacao_do_solo__cos_/cos__2007/
• ESRI. ArcGIS Desktop, Version 9.3. Redlands: Environmental Systems Research Institute, 2004.
• Fontes L, Tomé M, Coelho MB, Wright H, Luis JS, Savill P. Modeling dominant height growth of Douglas fir (Pseudotsuga menziesii (Mirb.) Franco) in Portugal. Forestry. 2004;76:509-23.
• Goovaerts P. Geostatistics for natural resources evaluation. Oxford: Oxford University Press; 1997. 483p.
• Goovaerts P. Geostatistical modeling of uncertainty in soil science. Geoderma. 2001;103:3-26.
• Godinho-Ferreira P, Azevedo A, Rego F. Carta da tipologia florestal de Portugal Continental. Silva Lusitana. 2005;13(1):1-34.
• Grant A, Worrell R, Wilson S, Ray D, Masonet B. Achieving diversity in Scotland’s forest landscapes. Forestry Commission Scotland Practice Guide. Forestry Commission, 2012. 38p.
• Hlásny T, Trombik J, Bosel’AM, Merganic J, Marusák R, Seben V et al. Climatic drivers of forest productivity in Central Europe. Agricultural and Forest Meteorology. 2017;234-35:258-73.
• Hock BK, Payn TW, Shirley JW. Using a geographical information system and geostatistics to estimate site index of Pinus radiate for Kaingaroa Forest. New Zealand Journal of Forestry Science. 1993;23:264-77.
• Huang SM, Titus SJ. An index of site productivity for uneven-aged or mixed-species stands. Canandian Journal of Forest Research. 1993;23(3):558-62.
• ICNF. IFN6 - Áreas dos usos do solo e das espécies florestais de Portugal continental. Resultados preliminares. Lisboa: Instituto da Conservação da Natureza e das Florestas; 2013. [acessado em: 05 jun. 2017] Disponível em: http://www.icnf.pt/portal/florestas/ifn/ifn6#dad
  » http://www.icnf.pt/portal/florestas/ifn/ifn6#dad
• ICNF. Cartografia nacional de áreas ardidas. Lisboa: Instituto da Conservação da Natureza e Florestas; 2017. [acessado em: 05 Jun. 2017]. Disponível em: http://www.icnf.pt/portal/florestas/dfci/inc/info-geo
  » http://www.icnf.pt/portal/florestas/dfci/inc/info-geo
• Kerr G. The use of silvicultural systems to enhance the biological diversity of plantation forests in Britain. Forestry. 1999;72(3):191-205.
• Lima JSS, Silva JTO, Oliveira RB, Almeida VS, Vanzo FL. Estudo da viabilidade de métodos geoestatísticos na mensuração da variabilidade espacial da dureza da madeira de paraju (Manilkara sp). Revista Árvore. 2006;30(4):651-7.
• Mandallaz D. Estimation of the spatial covariance in universal kriging: application to forest inventory. Environmental and Ecological Statistics. 2000;7:263-84.
• Meyer H.A. A Mathematical Expression for Height Curves. Journal of Forestry. 1940;38:415-520.
• Mitsuda Y, Satoshi I, Sakamoto S. Predicting the site index of sugi plantations from GIS-derived environmental factors in Miyazaki Prefecture. Journal of Forest Research. 2007;12:177-86.
• Nanos N, Calama R, Cañadas N, García C, Montero G. Spatial stochastic modeling of cone production from stone pine (Pinus pinea L.) stands in the Spanish northern plateau. In: Amaro A, Reed D, Soares P editors. Modelling forest systems. Wallingford: Cab International, 2003. p.131-41.
• Navalho I, Alegria C, Quinta-Nova L, Fernandez P. Integrated planning for landscape diversity enhancement,fire hazard mitigation and forest production regulation: A case study in central Portugal. Land Use Policy. 2017;61:398-412.
• Nunes L, Patrício M, Tomé J, Tomé M. Modeling dominant height growth of maritime pine in Portugal using GADA methodology with parameters depending on soil and climate variables. Annals of Forest Science. 2011;68(2):311-23.
• Oliveira A. Boas práticas florestais para o pinheiro bravo; manual. Centro Pinus; 1999. 32p.
• Oliveira A, Pereira JS, Correia AV. A silvicultura do pinheiro bravo. Centro Pinus; 2000. 111p.
• Páscoa F. Using forest inventory data to build growth and yield stands models. In: Wensel L, Biging G, editors. Proceedings of the IUFRO Conference on Forest Simulation Systems. University of California, Division of Agriculture and Natural Resources. 1990. p.279-86. (Bulletin 1927).
• Payn TW, Hill RB, Hock BK, Skinner MF, Thorn AJ, Rijkse WC. Potential for the use of GIS and spatial analysis techniques as tools for monitoring changes in forest productivity and nutrition, a New Zealand example. Forest Ecology and Management. 1999;122:187-96.
• Pelissari AL, Caldeira SF, Figueiredo Filho A, Amaral MS. Propostas de mapeamentos da capacidade produtiva de sítios florestais por meio de análises geoestatísticas. Scientia Forestalis. 2015;43(107):601-8.
• Pereira HG, Brito MG, Albuquerque T, Ribeiro J. Geostatistical estimation of a summary recovery index for marble quarries. In: Soares A editor. Geostatistics. 1993;2:1029-40.
• Pereira HG, Sousa AJ, Ribeiro JT, Salgueiro AR, Dowd P. Correspondence analysis as a modeling tool. 2015. 236p.
• Ribeiro A, Ferraz Filho AC, Tomé M, Scolforo JRS. Site quality curves for african mahogany plantations in Brazil. Cerne. 2016;22(4):439-48.
• Saito H, Goovaerts P. Geostatistical interpolation of positively skewed and censored data in a dioxin-contaminated site. Environmental Science and Technology. 2000;34(19):4228-35.
• Salas C. Site productivity of Nothofagus forests in Chile: a preliminary analysis. In: Lopes D, Tomé M, Liberato M, Soares P editors. Book of Abstracts IUFRO Conference 2010: Mixed and Pure Forests in a Changing World, 6-8 Oct. 2010, UTAD, 2010. p.119.
• Santos C, Almeida J. Caracterização espacial de um índice de produtividade nos povoamentos de pinheiro bravo em Portugal. Finisterra - Revista Portuguesa de Geografia. 2003;75:51-65.
• Skovsgaard J, Vanclay J. Forest site productivity: a review of the evolution of dendrometric concepts for even-aged stands. Forestry. 2008;81:13-31.

Publication Dates

• Publication in this collection
  2017

History

• Received
  07 June 2016
• Accepted
  24 May 2017