Developing an index for forest productivity... DEVELOPING AN INDEX FOR FOREST PRODUCTIVITY MAPPING – A CASE STUDY FOR MARITIME PINE PRODUCTION REGULATION IN PORTUGAL

1 Received on 07.06.2016accepted for publication on 24.05.2017. 2 Instituto Politécnico de Castelo Branco,Escola Superior Agrária,IPCB/ESA, Portugal. E-mail:<susanamestre@aflobei.pt>. 3 Centro de Estudos de Recursos Naturais, Ambiente e Sociedade,Instituto Politécnico de Castelo Branco,Escola Superior Agrária,CERNAS/IPCB/ESA, Portugal. E-mail: <crisalegria@ipcb.pt>. 4 Instituto Politécnico de Castelo Branco,Escola Superior de Tecnologia,IPCB/EST, Portugal. E-mail: <teresal@ipcb.pt>. 5 BioMedware, United States of America. E-mail: <Pierre.Goovaerts@biomedware.com>. *Corresponding author.

ALBUQUERQUE MTD et al.

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 3 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 (Alves et 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.

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 Developing an index for forest productivity... 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 2 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.
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): Eq1 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 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).

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  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 2 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 twostep 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(Pereira et al., , 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 Developing an index for forest productivity... 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.

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) backtransform 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).

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 50year 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 3 ha -1 yr -1 in low quality sites (Sh25<15), 8 m 3 ha -1 yr -1 in intermediate quality sites (Sh25Î [15,18[ m) and 12 m 3 ha -1 yr -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 ALBUQUERQUE MTD et al.
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.
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.

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 crossvalidation 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).

Site
Area of MAI 45 Relative Area of equal Long term Productivity Hectares to Years of productivity forest (m 3 ha -1 yr -1 ) productivity productivity sustainability equivalent* harvest per harvesting  (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.(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 "MAI 45 " 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.
Developing an index for forest productivity...
The PI map (Fig. 2a) after being reclassified in three productivity class (low -PI<-0.48;intermediate -PI Î ) 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 50year 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.

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 45 ) 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.

Legend: ( 1 )
Site productivity PI class -is classified as low for PI<-0.48;intermediate for PI Î [-0.48,-0.16[andhigh for PIe"-0.16;(2) Area of forest (ha) -is the maritime pine area by PI class evaluated using the PI map; (3) MAI 45 -is the mean annual increment of stand volume at the rotation age of 45 years (m 3 ha -1 yr -1 ) assigned to each PI class; (4) Relative productivity -is obtained as the ratio between "MAI 45 " in each PI class and the intermediate value of 8 m 3 ha -1 yr -1 ;

Figure 2 -
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).

Table 1 -
Area control method