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

Abstract

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 significantly influencing the productivity levels. 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 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.

1. 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].

Forest productivity estimations are necessary for both effective forest management decisions and the assessment of site conditions for ecological studies [10,11]. This knowledge greatly improves the ability to manage forest resources at a specific site, which in turn will help in making cost-effective decisions on the use of silvicultural soil [12]. Site productivity can be defined in many ways, depending on the objectives of resource managers [13]. From a forest management perspective, productivity is generally defined as the ability of a site to produce wood volume [14]. Assessing forest productivity is essential to understand sustainability in forest management because timber extraction should be based on the capacity for recovery of a forest area [15].

Commonly, forest monitoring has been carried out by traditional methods, which include extensive field sampling [16]. However, these methods are highly costly and time-consuming, especially in inaccessible locations [17,18]. Productivity assessments have traditionally been performed by determining the quality of the forest stands based on the site index (SI) approach [1]. For more than 100 years, the SI has been implemented as the world’s most widely used forest productivity measure. The SI based decision making serves as the baseline for land use management, forestry investment analysis, and growth and yield predictions [19,20]. The SI models are generally developed to describe the response of a defined variable, commonly the height of dominant trees [20,21], with respect to the quality of several edaphic, topographic, climatic and biotic variables [22].

For example, Ercanli et al. [12] developed a multiple linear regression model in Turkey for Picea orientalis with field sampling data and topographic, soil nutrients and climatic factors. The model explained 77% of the variation between factors and field data. In Ireland, Farrely et al. [23] generated a series of models to predict the SI of Picea sitchensis using climate variables and soil chemical factors as predictor variables. Bueis et al. [24] developed a discriminant model for estimating SI for Pinus yvestris L. in northern Spain using soil, climatic, and physiographic parameters. Sharma et al. [25] also developed a set of multiple linear regressions with different combinations of sites and climatic factors for Picea abies and Pinus silvestris. These studies agree that soil, topographic and climatic variables have the greatest influence on tree height and thus in SI modeling. In addition, studies have highlighted the feasibility of using accessible and low cost tools such as digital elevation models, lidar, or spectral images, to classify forest SI with high accuracy [10,11,26].

The use of SI curves based on height analysis is a preferred estimation of productivity because of the strong relationship between height and age of dominant trees [27], i.e., the age of healthy trees can serve as a predictor for the height of the trees. This and other plant-based analysis are called phytocentric methods while the geocentric methods are based on site properties such as climate, topography and soil [20]. It has been also stated that the productivity analysis should neither be constrained to the use of tree structural variables nor to the use of biophysical variables exclusively [20]. In many cases the use of both phytocentric and geocentric methods could be highly useful. Other procedures to estimate SI involve the use of some relationships such as the ratio of height/diameter as an inverse measure of productivity [28] or the use of non-parametric random forest models to estimate SI as a function of climatic variables [29]. In many cases, the SI cannot be measured directly and must be estimated from other site characteristics related to climate, topography, and soil [30]. However, the use of geographic variability is essential to account for local plant growth responses.

The increasing availability of spatially explicit data created by digital mapping using the geographic information system (GIS) allows the development of models with a spatial context for forest management that promise high regional accuracy. The combination of this spatial information with field data provides valuable information on the relationships between site condition and environmental variables [31]. For this reason, the joint GIS use of the SI is an attractive alternative in forest management [32,33]. The main benefit of the GIS is the capacity to model environmental problems, including large volumes of spatial data that link processes occurring across space and time [34]. Prediction models are applied to field data and as a result only provide point estimates of site productivity. The application of geostatistical techniques to layers of point data allows to pass 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.

2. Materials and Methods

2.1. 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. Summary of site and forest characteristics of the P. arizonica sampled plots of the study.

Variable	Average	Max Value	Min Value
Above sea level altitude (m)	2314	3020	1690
Slope (%)	31	85	5
A-Horizon (cm)	33	97	3
Rockiness (% Cover)	23	78	0
Bare soil (% Cover)	13	85	0
Litter depth (cm)	2	9.5	0.3
Min temperature °C	5	13	1
Max temperature °C	23	32	17
Precipitation (mm)	898	1248	576

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

2.2. 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 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² 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² 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 (1940–2004). 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.

2.3. Site Index Model

In this study, the Schumacher model (Equation (1)) (see Palahi et al. [50] for details) was evaluated.

$$H={\beta }_0{e}^{-{\beta }_{1}E-1}$$

(1)

where H is the dominant height (m); E is the age (years); and β₀ and β₁ 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², and the distribution of residuals.

2.4. 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]

$$K_{\mathit{APPA}}=\frac{N{{\displaystyle \sum }}^K{X}_{ii}{{\displaystyle \sum }}^k\left(X_{i+}\ast X_{+i}\right)}{N^2-{{\displaystyle \sum }}^k\left(X_{i+}\ast X_{+i}\right)}$$

(2)

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]. These authors classified the effectiveness as poor (0–0.4), fair (0.4–0.55), good (0.55–0.7), very good (0.7–0.85), and excellent (0.85–1.0).

3. 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. Descriptive statistics of the sampled trees.

	DBH (cm)				Height (m)
Age	Mean	St. Dev.	Min	Max	Mean	St. Dev.	Min	Max
10	4.3	2.5	3.1	10.6	2.5	0.9	1.6	4.6
20	8.9	4.2	5.5	16.6	5.5	1.9	2.1	9.7
30	13.8	5.1	4.3	21.6	8.8	2.8	3.5	12.5
40	18.7	5.3	9.3	26.7	11.7	3.0	5.7	16.6
50	22.6	5.4	12.6	32.2	13.6	2.9	7.2	18.7
60	26.3	5.1	16.0	34.4	15.5	3.0	9.6	19.8
70	28.1	5.3	20.1	35.0	16.4	3.5	11.3	23.9
80	31.4	6.4	23.9	40.9	17.8	3.4	12.1	25.0
90	33.1	6.1	26.2	43.7	19.6	3.5	14.7	25.6

DBH = Diameter at Breast Height, St.Dev. = Standard Deviation, Min = Minimum, Max = Maximum.

. 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. Statistical fit for the Schumacher model, parameters values and their significance.

Schumacher	SSE	MSE	RMSE	R2	R2ajus	DW	Parameter	Estimate	Approx Std Err	t Value	Approx Pr > |t|
	554.6	2.6287	1.6213	0.9258	0.9247	1.3877	β0	25.5103	1.0893	23.42	<0.0001
							β1	31.7767	1.8258	17.40	<0.0001
							p1d	0.9656	0.0208	46.39	<0.0001
							p2d	0.8785	0.0251	34.95	<0.0001

SSE = Sum Square Error, MES = Mean Square Error, RMSE = Root Mean Square Error, DW = Durbin Watson Statistic, Approx. Std Err = Approximate Standard Error, Approx. Pr > |t| = Approximate p value.

). 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. Mean values of physical variables in areas with different productivity levels of P. arizonica in the area of Guadalupe y Calvo, Chihuahua, Mexico.

Variable	Statistics				Productivity Level
					Low		Medium		High
	SS	MSE	Fc	Pr > F	Mean	Stdv	Mean	Stdv	Mean	Stdv
Slope (degree) **	0.215	0.107	3.999	0.020	37A	15	30A	15	12B	8
Soil depth (cm) **	527.354	263.677	0.698	0.012	26A	18	36B	19	44B	11
Rockiness (%)	0.154	0.077	2.332	0.100	29	2	22	17	21	18
Bare soil (%)	0.002	0.001	0.056	0.945	12	1	13	14	13	17
A-Horizon (cm)	4.311	2.156	2.039	0.133	1.5	0.8	1.6	0.8	1.9	1.3
Litter depth (cm)	0.071	0.035	0.644	0.527	0.6	0.2	0.6	0.2	0.6	0.3
Minimum temperature (°C) **	37.378	18.689	3.261	0.041	5A	2	4A	1	8B	3
Maximum temperature (°C) **	86.152	43.076	5.208	0.006	23A	3	23A	2	18B	3
Precipitation (mm) *	1219.023	609.511	0.043	0.096	794A	43	897B	92	901B	116

SS = Sum of Squares, MSE = Mean Square Error, Fc = Calculate F, Pr > F = p value, Stdv = Standard Deviation, * Variables that were significantly different among productivity levels at p ≤ 0.1, ** Variables that were significantly different among productivity levels at p ≤ 0.05.

).

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.

4. 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² = 0.93) our results are comparable to SI models obtained for a hybrid Larix x eurolepis in Sweden (R² = 0.99) [63], P. silvestris in Spain (R² > 0.92) [50], and other pine species in Mexico (R² > 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² 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).

5. 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 diversity, forest exploitation or anthropogenic activities like agriculture and livestock, on the productive variability of forest species.