A Spatial Forestry Productivity Potential Model for Pinus arizonica Engelm, a Key Timber Species from Northwest Mexico

: Pinus arizonica is a widely distributed tree species growing in temperate forests of Northwest Mexico where it is utilized through different regeneration harvest methods. Yet, management models based on estimations of its productive potential are sorely lacking. In this study, a procedure to create a productive map using site index (SI) equations and Geographic Information Systems (GIS) was developed. A SI model for P. arizonica was created for the study area and used to classify a group of randomly sampled plots on three productivity categories (High, Medium, and Low) for management purposes. Climatic, topographic and edaphic variables were determined on the sampled plots. Then, a statistically-based analysis was performed to identify the climatic, topographic and edaphic variables signiﬁcantly inﬂuencing the productivity levels. Based on the values of these signiﬁcant variables, a map of productive potential was elaborated for the whole study area. Sites with the highest productivity were those with slopes ≤ 12 ◦ , soil depths ≥ 0.46 m, minimum and maximum mean annual temperatures of 5 ◦ C and 18 ◦ C respectively, and precipitation ≥ 900 mm. This methodology could be considered for similar species/conditions where productivity models do not exist or to update old models rendered obsolete by climate change.


Introduction
Forest ecosystems play an important role in biodiversity conservation and the provision of ecosystem services, such as timber production and carbon sequestration [1]. These ecosystems play an important role in the economies of worldwide rural communities [2,3]. The need to provide forest products to a growing world population despite major obstacles such as climate change, desertification, environmental pollution, and loss of biodiversity is a major challenge for forest managers [4]. Forest ecosystems have experienced strong degradation due to deforestation, changes in land use and land cover, fires, climate variability and fragmentation [5][6][7][8]. Forests in northern Mexico are used extensively for timber extraction but have experienced disturbance by changes in land use by agricultural and pastoral activities, resulting in productivity changes as well as habitat modification for many species of fauna and flora [9]. from a point representation to a continuous one, as well as the obtaining of the uncertainty measure of the sites not sampled [35]. Several studies have been developed to create useful maps for forest management [36][37][38][39]. However, attempts to map geospatial attributes to measure forest productivity have been limited [37]. To bridge the gap between geospatial information, the relationship with forest productivity and reduce fieldwork, SI and GIS as a set of spatial variables such as climate, soil, and topography, allow the development of maps, which could be effective tools for managing complex ecosystems [40][41][42]. Recent approaches to bridge the gap between geospatial information and its relationship to SI have been employed by Waring et al., (2014) [40], where they used a forest growth model based on 3-PG (Physiological Principles Predicting Growth) processes using satellite estimates of the Maximum Foliar Area Index contrasting it with the site's growth potential for Pseudotsuga menziesii and thus visualizing its variation in the western region of North America. Brandl et al., (2014) [31] used national forest inventory information from the state of Bavaria, Germany, together with maps derived from the digital elevation model, temperature and precipitation maps, soil maps, to spatially predict the productivity of the SI, using the Generalized Additive Model. Mestre et al., (2017) [33] used geospatial slope and exposure variables together from a productivity index map constructed using a multi-Gaussian Kriging method in a geographic information system, where they evaluated high, intermediate, and low productivity levels.
Despite the development of allometric equations and fitting models for SI [43], the incorporation of geospatial variables influencing tree growth and the use of GIS tools have been poorly studied and applied to forests management in Mexico. This is particularly the case of the temperate forests of Northwest Mexico, where Pinus arizonica Engelm is a species distributed on a wide range of soil and topographic conditions [44]. The need for current forest productivity estimations, not only in Mexico but worldwide, arises also because of contemporary climate change, as forest growth estimations based on decades-old equations may not reflect current tree growth rates [40].
Pinus arizonica is a key species of great ecological, economic, and social importance in the Sierra Madre Occidental, which includes the states of Chihuahua and Durango, Mexico [45,46]. Yet, little information exists on its productivity related to the spatial environmental variability and thus timber utilization cannot have a sustainability target. Here, we attempted to analyze the validity and usefulness of growth curves linked to GIS techniques to generate an easy-to-use geographic model to predict tree productivity. The objectives of this study were (1) to develop a SI model for P. arizonica under the conditions of southern Chihuahua, Mexico, considering the height and age of dominant trees, and (2) to create a spatial model to display the landscape productivity for this species as a function of topographic, climatic and edaphic variables on a scale relevant to regional forest management. We used the Schumacher model as basis for our SI analysis and the obtained productivity classes were assigned according to the conditions of physical variables in the field to obtain useful predictor maps. The results could be used by forest managers and planners to better understand forest productivity and to implement decisions based on sustainable forest management plans. Likewise, geospatial information could be useful for territorial ecological ordering considering the spatial distribution of forest productivity.

Study Area
This study was conducted in the area of Guadalupe y Calvo, Chihuahua, which is part of the Sierra Madre Occidental, in Northwest Mexico ( Figure 1). This area was chosen because P. arizonica is a key species naturally and amply distributed in these forest ecosystems. The altitude of the study area ranges from 1600 m to 3020 m, with an average slope of 33% (Table 1). The most common soil types are regosols and lithosols with an average A-horizon depth of 33 cm, about 23% of rockiness, and an average mulch depth of 2 cm. The predominant climate is temperate sub-humid with an annual precipitation between 576 mm to 1248 mm, mainly occurring as rainfall in the summer, and with occasional snowfall occurring from December to February. In this region, the main productive activities are timber harvesting, cattle ranching, rainfed agriculture, and mining. The forests of this region are of great economic importance and, particularly, P. arizonica is one of the most exploited species due to its timber quality [47].

Methodological Outline and Data
To create a spatial model of productivity we used physical variables of soil, topography and climate. The general methodological outline to elaborate such spatial model appears in Figure 2. The information was obtained from pure and mature stands of P. arizonica distributed over different aspects. We selected 42 P. arizonica dominant trees (tallest trees of the stands) from different conditions of slope, soil depth and aspect in the study area. Then, we carried out a stem analysis following the methodology proposed by the British Columbia Forest Productivity Council [48]. Each

Methodological Outline and Data
To create a spatial model of productivity we used physical variables of soil, topography and climate. The general methodological outline to elaborate such spatial model appears in Figure 2. The information was obtained from pure and mature stands of P. arizonica distributed over different aspects. We selected 42 P. arizonica dominant trees (tallest trees of the stands) from different conditions of slope, soil depth and aspect in the study area. Then, we carried out a stem analysis following the methodology proposed by the British Columbia Forest Productivity Council [48]. Each tree was harvested to subsequently perform cross-sectional cuts at 0.30 m, 1.30 m (diameter at breast height or DBH) and every 2 m until the stem reached 10 cm in diameter. The length of the remaining tip was recorded. Following guidelines, we counted the number of annual rings, measured diameters of the tree was harvested to subsequently perform cross-sectional cuts at 0.30 m, 1.30 m (diameter at breast height or DBH) and every 2 m until the stem reached 10 cm in diameter. The length of the remaining tip was recorded. Following guidelines, we counted the number of annual rings, measured diameters of the stem transverse cross sections, and estimated the height of the tree at different ages. This information integrated the database of height-age data pairs to perform the SI analysis. To identify the influence of the physical variables over the forest productivity 220 circular plots of 1000 m 2 each with the presence of P. arizonica were located and sampled by using a 5 km systematic design. Every site was georeferenced using a Garmin Etrex Global Positioning System (GPS) and then the height and the age of the two dominant trees were registered for each plot. The height was measured with a Suunto clinometer, whereas the age was determined by counting the trees annual rings through samples extracted with a Pressler drill. At the center of each sampling plot we recorded slope and aspect and located a 1 m 2 sub-plot to determine percent bare soil area, pedregosity, and litter depth. Also, a soil pit was dug adjacent to the subplot to determine soil depth. We also included in the database climate variables such as the means of minimum and maximum annual temperatures and annual precipitation from a time series of 64 years . Climate data were obtained from ERIC III, a historical weather information depository provided by the Mexican Meteorological Service (Mexican Institute of Water Technology) [49] (https://www.imta.gob.mx/productos/software/eric-iii-version-3-2-extractor-rapido-deinformacion-climatolo-detail). The selected climate data corresponded to weather stations placed around and within the study area. To identify the influence of the physical variables over the forest productivity 220 circular plots of 1000 m 2 each with the presence of P. arizonica were located and sampled by using a 5 km systematic design. Every site was georeferenced using a Garmin Etrex Global Positioning System (GPS) and then the height and the age of the two dominant trees were registered for each plot. The height was measured with a Suunto clinometer, whereas the age was determined by counting the trees annual rings through samples extracted with a Pressler drill. At the center of each sampling plot we recorded slope and aspect and located a 1 m 2 sub-plot to determine percent bare soil area, pedregosity, and litter depth. Also, a soil pit was dug adjacent to the subplot to determine soil depth. We also included in the database climate variables such as the means of minimum and maximum annual temperatures and annual precipitation from a time series of 64 years . Climate data were obtained from ERIC III, a historical weather information depository provided by the Mexican Meteorological Service (Mexican Institute of Water Technology) [49] (https://www.imta.gob.mx/productos/software/eric-iii-version-3-2-extractor-rapido-de-informacion-climatolo-detail). The selected climate data corresponded to weather stations placed around and within the study area.

Site Index Model
In this study, the Schumacher model (Equation (1)) (see Palahi et al. [50] for details) was evaluated.
where H is the dominant height (m); E is the age (years); and β 0 and β 1 are empirical parameters of the model. Our analysis considered a tree base age of 80 years. This model was tested using the guide curve (GC) methodology [27]. To evaluate the results of the model we considered the value of the root mean square error (RMSE), the value of R 2 , and the distribution of residuals.

Analysis
The analysis to determine the best-fitted model was performed through the NLIN procedure and the DUD method by using the Statistical Analysis System version 10.1 (SAS Inc ® Cary, NC, USA). The Shapiro-Wilk test was conducted to verify compliance with the assumptions of regression. Additionally, the relative and cumulative residual frequencies were obtained. These frequencies resemble a straight line with respect to the probability of normal distribution and the corresponding percentages resemble a Gauss-bell curve [51].
We used the predicted height according to the model using the guide curve method, and then a group of anamorphic site index curves were fitted. Tree growth was stratified in three quality classes (25 m = High, 21 m = Medium and 12 m = Low). The curves showing the growth trend were constructed by holding the shape parameters of the chosen model constant and varying the asymptote parameter as necessary to achieve the required dominant height when the tree age equaled base age, which was 80 years in this study. This approach has been used in previous studies [52].
The three height classes generated with the Schumacher model were used to assign the quality level of productivity to each of the 220 plots according to the height-age data of their dominant trees. Then, we performed analyses of variance (ANOVA) to determine if each of the physical variables recorded from the plots significantly influence the productivity level (High, Medium, and Low). A comparison of means with the Tukey test (α = 0.1) was also performed to determine differences of such variables among productivity levels.
In a following step, a spatial layer for each of the physical variables significantly influencing the productivity level was created with the software ArcGIS 10.2 [53]. The layer of topography (slope) was elaborated from a Digital Elevation Model (DEM) 1:50,000 [54]. The layer of soil depth was taken from a soil digital map of the state of Chihuahua, Mexico. The climatic layers (precipitation and temperatures) were created by interpolating the data of climate variables obtained from ERIC III. Interpolation was based on the Inverse Distance Weighting method [55], employing the Geostatistical Wizard available in the software ArcGIS 10.2 [53]. The cell size used for all the layers was of 25 m.
For each layer, their values were categorized in three levels of productivity; High (3), Medium (2), or Low (1), according to the results from the Tukey tests. Then, the mean value of each level was assigned to the pixels based on their corresponding category. Once categorized, the layers were summarized (slope + soil depth + mean minimum temperature + mean maximum temperature + precipitation). As a result, a map was created, which was then reclassified in three equal intervals (High, Medium, and Low) to create the productive potential map.
To validate this new map model, age and height of the trees were taken from a random sample of 190 sites within the study area. The effectiveness of the productive potential map was assessed with the K APPA Index, which is denoted by Equation (2) [56] where: K APPA = Kappa index, k = number of matrix files, x ii = observation number of row i and column i (along the diagonal), x i+ and x +i = Total marginal for row i and column i, respectively, and N = total number of observations. To characterize the level of agreement suggested by the K APPA index, we used the qualitative descriptors pointed out by Monserud and Leemans [57].

Results
Descriptive statistics of the diameter at breast height (DBH) and height at different ages of the dominant trees used to fit the SI model are shown in Table 2. We considered records of every ten years. Such records were gotten from a stem analysis performed by following the methodology proposed by the British Columbia Forest Productivity Council [48]. The maximum height registered from the trees was 25.6 m at an age of 90 years, with diameters of 33.1 ± 6.1 cm, which may indicate the existence of high variability in the productivity of the study area. This may be an effect of the high topographic and soil variability of the sites where the trees grew. This response could be observed in young (10-30 years), as well as in full mature trees (50-90 years). The Schumacher model [50] showed an adequate fit to the stem analysis data (Table 3). Although we initially found positive autocorrelation in the residuals, estimated with the Durbin-Watson statistic (DW = 0.33), the autocorrelation was corrected (DW = 1.4) by applying two delays in the residuals using a CAR (2) shape structure [58]. The quality of the SI model adjustment for P. arizonica was then in compliance with the assumptions of normality of the regression sample errors, since the SW test was acceptable with W normal = 0.99. In addition, no trend was observed on the variance of the residual distribution, thus discarding heteroscedasticity. The resulting equation used to classify the SI on the study area was H = 25.21029e −31.7766 E−1 . By using this equation, the guide curve was estimated considering a base age of 80 years. Upon defining the base age and the trend of the guide curve, the SI was established at 25 m, 21 m, and 12 m for High, Medium, and Low levels of productivity, respectively. This fitted model produced the height curves shown in Figure 3. These curves were plotted together with the observed data. The shape of the curves is very close to the shape of the observed data with realistic asymptotes and growth patterns. This allowed us to assign the levels of productivity to the group of plots sampled on the field.
The ANOVA results showed statistical differences among productivity levels for slope, precipitation, minimum and maximum temperature, and soil depth. As previously stated, autocorrelation was not a limiting factor for our result analysis. The values associated to the highest productivity level of P. arizonica through the selected model were 12 • of slope, 46 cm of soil depth, 901 mm of precipitation, as well as 5 • C and 18 • C of minimum and maximum mean annual temperature, respectively (Table 4).
Following the methodology, the map of productivity was created, and an estimation of the surface for each productivity class was made. The sites with High productivity in the area occupied by the temperate forest corresponded to 6.5% (43,112 ha), whereas the sites with Medium and Low productivity corresponded to 60.4% (402,887 ha), and 33.1% (220,679 ha) of the study area, respectively ( Figure 4). The rest of the surface (26%) (236,480 ha) is considered as non-forest area and corresponded to deciduous forest with an elevation lower than 1800 m, where the presence of P. arizonica was not detected. The accuracy of the spatial model was assessed by comparing values from the map with records of productivity measured on the field. The estimated value for K APPA Index was 0.79 ± 0.08, resulting in a very good agreement.  Figure 3. These curves were plotted together with the observed data. The shape of the curves is very close to the shape of the observed data with realistic asymptotes and growth patterns. This allowed us to assign the levels of productivity to the group of plots sampled on the field. The ANOVA results showed statistical differences among productivity levels for slope, precipitation, minimum and maximum temperature, and soil depth. As previously stated,

Discussion
Estimating the distribution of the productive potential through traditional sampling is a costly and time-consuming activity [59]. In this study we show a parsimonious methodology that allows mapping the spatial distribution of the productive potential of P. arizonica which has been carried out in other countries such as Germany [31], Spain [32], and Finland [38]. Forest productivity estimations are essential to achieve a sustainable forest management [60]. The most utilized indirect method worldwide is modeling height of dominant trees [61]. This is generally performed with data pairs of age-height from stem analysis of trees growing in permanent monitoring plots [62]. The use of data from stem analysis has been extensively documented [61]. Thus, we used stem analysis data to fit the P. arizonica SI model. Our study constitutes a methodological proposal for mapping and documenting the SI in the northwestern region of the state of Chihuahua, using field information and GIS. Because of the high adjustment value that we obtained, our ability to map local information (weather stations) and available spatial information (digital elevation model, soil map), the proposed approach confirms the successful use of GIS, field data sets, and overlapping of layers. This offers the possibility for land managers to obtain and generate accurate information at the scales of their interest and with information from their inventories.
Our results showed an adequate growth prediction for all tree ages, which ratifies the effectiveness of this procedure in our study area. Since acceptable statistical and graphical results

Discussion
Estimating the distribution of the productive potential through traditional sampling is a costly and time-consuming activity [59]. In this study we show a parsimonious methodology that allows mapping the spatial distribution of the productive potential of P. arizonica which has been carried out in other countries such as Germany [31], Spain [32], and Finland [38]. Forest productivity estimations are essential to achieve a sustainable forest management [60]. The most utilized indirect method worldwide is modeling height of dominant trees [61]. This is generally performed with data pairs of age-height from stem analysis of trees growing in permanent monitoring plots [62]. The use of data from stem analysis has been extensively documented [61]. Thus, we used stem analysis data to fit the P. arizonica SI model. Our study constitutes a methodological proposal for mapping and documenting the SI in the northwestern region of the state of Chihuahua, using field information and GIS. Because of the high adjustment value that we obtained, our ability to map local information (weather stations) and available spatial information (digital elevation model, soil map), the proposed approach confirms the successful use of GIS, field data sets, and overlapping of layers. This offers the possibility for land managers to obtain and generate accurate information at the scales of their interest and with information from their inventories.
Our results showed an adequate growth prediction for all tree ages, which ratifies the effectiveness of this procedure in our study area. Since acceptable statistical and graphical results were obtained with the Schumacher model in this study, testing more equations was considered not necessary. Due to the close relationship between tree height and tree age that we found for P. arizonica (R 2 = 0.93) our results are comparable to SI models obtained for a hybrid Larix x eurolepis in Sweden (R 2 = 0.99) [63], P. silvestris in Spain (R 2 > 0.92) [50], and other pine species in Mexico (R 2 > 0.96) [43]. The statistical analysis clearly showed that the SI is adequate to measure forest productivity as indicated by the study carried out by Günlü et al. [64] in an unmanaged forest where analyze the relationship between direct and indirect methods using the Chi-square test. The test indicated a statistically significant relationship between the SI determined by the indirect method is satellite imagery.
Even though high levels of agreement between the real and predicted values can be achieved with the SI models, they still present some limitations. Such limitations include that these models are developed and are effective for monospecific and even age stands. In addition, for the development of these models it is supposed that the quality of the stand is calculated from the information of the specific points in which the sampled trees are located; however, it does not include the variability of the physical or biological variables all over the stand [65] and geographical location between ecoregions [32].
Moreover, thinning or hard grading tend to modify the height of the trees growing in the stand, which also modifies the dominant height causing over or sub estimations on the trend of the trees growth. In these cases, the use of soil-site relationships is more suitable for estimating site quality [23]. With these limitations in mind our modeling results had adequate agreement with our ground-truth comparisons [57]. The visual interpretation of the SI and the map of productive potential revealed an expected agreement with the topographic variations and the distribution of the SI, found in the study area. The High SI level appeared near urban-rural areas, indicating the presence of management in the forests, while the Low productive potential SI appear in the most remote areas, where the forests are possibly unmanaged, affecting the quality of the SI [66].
Our results highlight the importance of relating SI models with direct methods through the use of biophysical variables to generate continuous models that are acceptable and easy to understand. We integrated a procedure combining both direct and indirect methods to generate spatial maps of productivity. We accomplished that by classifying the biophysical variables influencing the height of dominant trees. Our findings showed that slope, precipitation, minimum and maximum temperature, as well as soil depth are the five variables with the highest relationship with the SI model of P. arizonica. The spatial distribution of the SI was consistent with the gradient of temperature and precipitation in the municipality. Therefore, an increase in SI productivity is probably the result of increased precipitation and low temperatures which is documented by Reich et al. [67] and Peters et al. [68].
The approach of this study shows that field work together with digital maps of slope, precipitation, minimum and maximum temperature, as well as soil depth provide a representative sample of the relationship between average forest productivity across the SI and the variability of biophysical conditions as reported by Waring et al. [40] at least on a regional scale [69].
That allowed elaborating a productivity map with a good accuracy. Productivity maps as developed with our procedure can be a useful and handy tool for forest managers as the maps show spatial patterns easy to visualize. Besides, maps created by these procedures could allow to estimate forest productivity even in areas where no forest inventories exist [68].
The values associated to the highest productivity level of P. arizonica, through the selected model were 12 • of slope, 46 cm of soil depth, and 5 • C and 18 • C of minimum and maximum annual temperature, respectively. Our results agree with other findings obtained for forest areas of Mexico, where productivity is commonly affected by slope, aspect and soil depth [70]. In agreement with this, it has been reported that soil variables may explain about 55% of the variation in forest productivity [71].
Other studies have found influence by other more specific variables such as soil pH and organic matter [61]. This is due to the scale of work we use, which is a regional scale [72]. If more detailed investigations are conducted, probably some other micro variable could be found influencing the growth of the trees. Yet, when trying to incorporate micro variables on mapping procedures, map generation could be prohibitively expensive and time consuming [71].
Multiple regression techniques and multivariate analysis have also been effective in the evaluation of forest productivity when employing direct methods because they allow inferring on the influence of complex variables over the variability of the productivity phenomenon [73]. Some examples of these are studies developed for Picea in forest of Turkey [13], which reported R 2 values of 0.77 when variables of soil, climate and topography were included in the analysis. Moreover, another study carried out also in Turkey compared three methods for the estimation of productivity: (1) indirect model of dominant height with guide curve, (2) direct, using soil samples, and (3) using satellite images. Differences among the methods were non-significant and the authors recommended the use of direct methods, especially in degraded and open land areas [10].
Our findings confirm the importance of using direct or combined methods to evaluate productivity of dry forest, such as the forest in northern Mexico. The use of maps of productivity could also offer reliable and easy to apply tools for the local forest managers. Additionally, the maps of productivity, which were created by using the influence of biophysical variables over species growth, could also allow to spatially classifying the productivity of lands without presence of trees. That could be useful for multiple purposes such as the establishment of forest plantations with the studied species [74,75]. However, the establishment of monitoring plots could be valuable to evaluate the productive variability of the species with a greater precision, also considering some factors such as the genetic variability [64] or the influence of forest, agricultural or livestock exploitation activities. Overall, we agree with Rodhouse et al. [76] that biogeographic modeling is a parsimonious approach to connect land management decisions with spatial monitoring.
Based on our results, we recommend the use of variables such as slope, temperature, and soil depth to perform cartographic studies and to estimate the productivity of these types of ecosystems in similar parts of Mexico and the world. The use of outdated growth forest models may be risky considering that global modification of climate is altering current vegetation growth patterns [40]. We expect our results and models to provide sorely needed tools that can be adopted for forest management in other regions of Mexico as part of the national strategy to increase forest productivity and production (Estrategia Nacional de Manejo Forestal Sustentable para el Incremento de la Producción y Productividad-ENAIPROS).

Conclusions and Recommendations
The combination of phytometric (SI models) and geocentric analysis provides a powerful tool for forest productivity. The use of variables such as slope, temperature, and soil depth to perform cartographic studies aimed at estimating forest productivity in temperate, mountainous, sub-humid areas was supported by our results. Based on the values of these significant variables, a map of productive potential was elaborated for the whole study area. Sites with the highest productivity have slopes ≤12 • with soil depths ≥0.46 m, with minimum and maximum mean annual temperatures of 5 • C and 18 • C respectively, and precipitation ≥900 mm (α ≤ 0.1).
We recommend to generate growth forests models where they do not exist and to update old models following the methodology we proposed here. The use of outdated growth forest models involves the risk of assuming growth patterns that do exist any longer because of environmental changes imposed by climate change. We expect our results and models to provide sorely needed tools that can be adopted for forest management in Mexico.
Additionally, these results could be used to identify sites for the potential establishment of permanent monitoring plots to learn more about the dynamics of forest growth. For further research, we recommend to carry out more detailed studies about the influence of variables such as genetic