Patterns of Density and Production in the Community Forests of the Sierra Madre Occidental, Mexico

: The Mexican Sierra Madre Occidental (SMO) represents a region where hundreds of plant species reach the limits of their northern or southern range. The SMO also features a unique cultural diversity, and many communities living within the forest or in its close vicinity depend on the products and services that these forests provide. Our study was based on a large set of remeasured ﬁeld plots placed in the forests of Durango which are part of the SMO. Using hierarchical clustering, three distinctly di ﬀ erent forest types were identiﬁed based on structural di ﬀ erences and the relation between stem density and basal area. Maximum forest densities were estimated using a 0.975 th quantile regression. Forest production (expressed as current periodic volume increment per unit of area and time) was estimated based on number of stems, forest density, mean height, and forest diversity. Forest density is the principal factors a ﬀ ecting periodic volume production. The discussion presented recommendations for the sustainable use of this unique natural resource. Maintaining minimum levels of residual density is key to ensuring the continued viability of the forests of the Mexican SMO. Future research is needed to identify optimum residual structures, productive residual densities, and desirable levels of biodiversity.


Introduction
Forest density is a multifaceted phenomenon, and the terminology relating to forest density is often the cause of misunderstanding. To prevent imprecise and confusing usage of the term, there is a need for clarification. In this study, we defined forest density as the degree of site occupancy by trees, i.e., the total tree biomass per unit forest area. It is almost impossible to assess that quantity. Therefore, basal area (the sum of the cross-sectional areas of the trees at breast height) is often used as a practical substitute for biomass. Forest density can be expressed in absolute or relative terms. Absolute measures of density are determined directly for a given location. Relative density is used to compare a given density with some standard, such as the full stocking of a yield table [1,2].

Permanent Field Plots
The observations used in this study were obtained from a network of permanent sample plots that are used to monitor the dynamics of the forests of Durango [19]. Altogether, 217 undisturbed plots without intermediate harvesting between remeasurements were selected for this study. a summary of the plots is presented in Table 1 for the two inventories. The plot establishment began in 2007, and remeasurements started in 2012. Each plot was remeasured after five years. The research plots included the main types of the forest stands in the SMO, providing an opportunity to study patterns of forest structure, density and production in this unique socioecological region. Each plot covered an area of 50 m × 50 m. The plots were distributed systematically, with a variable grid ranging from 3 km to 5 km depending on the size of the land properties. The tree number, species, dbh, total tree height (h; m), height to the live crown (hc; m), azimuth (A; degrees), and distance (d; m) from the center of the plot of all trees equal or greater than 7.5 cm in dbh were recorded.
The forests of the SMO have been managed by local communities for more than 100 years, mainly using selective harvest for sustainable timber production, but also for the maintenance of biological diversity and uneven-aged stand structures. The Continuous Cover Forestry System is preferred in the region because it promotes mixed and irregular forest stands and high biodiversity. It is known as the Mexican Method of Irregular Forest Regulation and is characterized by the absence of a rotation age that defines the time of the harvest. The stand age is undefined because trees of all ages usually occur in close vicinity to each other. Commercial harvests are based on maintaining growing stock levels within an ideal and negative exponential diameter class distribution, which is also known in the literature as an inverse J-shaped distribution. In uneven-aged stands, the idealized diameter distribution is thought to be balanced, meaning that it can be maintained by applying a given harvest rate in perpetuity [16].

Stratification of Field Plots
The distribution of diameters and heights were estimated by the Weibull model [20] using the function fitdistr of the MASS package [21] of R 3.6.0 [22]. The Weibull individual scale and shape plot parameters estimates were then grouped using the R function hclust [22] to obtain, iteratively and based on the method of centroids, a hierarchical cluster of the parameter estimates [23]. Among the different choices of cluster analysis, the methodology implemented in the hclust function of R, which is based on the complete linkage method for hierarchical clustering, was applied. This method defines the cluster distance between two clusters to be the maximum distance between their individual components. At every stage of the clustering process, the two nearest clusters were merged into a new cluster. The process was repeated until the whole data set was agglomerated into one single cluster.
In addition to careful inspection of the dendrogram classification of clusters, the Kruskal-Wallis Rank Sum Test (a nonparametric equivalent to ANOVA) was used to confirm or reject the results of the cluster process, specifically to test whether the sample plots came from unrelated population, or whether the population distributions were identical without assuming them to follow the normal distribution. The Kruskal-Wallis Rank Sum Test was applied to further discriminate among different groups of field plots, involving the plots for the variables of N, G, Dq, Ho, PAI, and S.
The plots were divided into three groups, which differed in terms of basal area, diameter distribution, and the relationship between number of trees and average tree size ( Figure 1). There was no obvious geographical preference. Each of the three groups of plots was widely distributed within the state of Durango, and the groups often mingled at close range.
Forests 2020, 11, x FOR PEER REVIEW 4 of 14 new cluster. The process was repeated until the whole data set was agglomerated into one single cluster.
In addition to careful inspection of the dendrogram classification of clusters, the Kruskal-Wallis Rank Sum Test (a nonparametric equivalent to ANOVA) was used to confirm or reject the results of the cluster process, specifically to test whether the sample plots came from unrelated population, or whether the population distributions were identical without assuming them to follow the normal distribution. The Kruskal-Wallis Rank Sum Test was applied to further discriminate among different groups of field plots, involving the plots for the variables of N, G, Dq, Ho, PAI, and S.
The plots were divided into three groups, which differed in terms of basal area, diameter distribution, and the relationship between number of trees and average tree size ( Figure 1). There was no obvious geographical preference. Each of the three groups of plots was widely distributed within the state of Durango, and the groups often mingled at close range. Figure 1. Spatial distributions of field plots by group in the state of Durango. The plots were divided into three groups, which differed in terms of basal area, diameter distribution, and the relationship between number of trees and average tree size. Plot affiliation to a particular group is indicated by one of three groups.

Maximum Forest Density
Populations of trees growing at high densities are subject to density-dependent mortality or selfthinning. For a given average tree size, there is a limit to the number of trees per area that may coexist in an even-aged stand [24]. A self-thinning model was used to describe the allometric relationship between the average tree size and the number of trees per hectare [4,[25][26][27]. The limiting relationship was described by the following exponential equation: where is the maximum number of living trees per ha (i.e., the full stocking density level), ln is the natural logarithm, Dq indicates the quadratic mean diameter (cm), α is an empirical coefficient, and β is the slope of the relationship between lnN and lnDq in fully stocked stands. The maximum density was estimated using the 0.975th quantile [28] of the R function quantreg [29]. In addition, we also calculated the Reineke's Stand Density Index (SDI) for each plot (using our own empirical exponent 1.536, instead of Reineke's 1.605). The SDI has been used as a measure of density for predicting forest productivity [30]: Figure 1. Spatial distributions of field plots by group in the state of Durango. The plots were divided into three groups, which differed in terms of basal area, diameter distribution, and the relationship between number of trees and average tree size. Plot affiliation to a particular group is indicated by one of three groups.

Maximum Forest Density
Populations of trees growing at high densities are subject to density-dependent mortality or self-thinning. For a given average tree size, there is a limit to the number of trees per area that may co-exist in an even-aged stand [24]. a self-thinning model was used to describe the allometric relationship between the average tree size and the number of trees per hectare [4,[25][26][27]. The limiting relationship was described by the following exponential equation: where N max is the maximum number of living trees per ha (i.e., the full stocking density level), ln is the natural logarithm, Dq indicates the quadratic mean diameter (cm), α is an empirical coefficient, and β is the slope of the relationship between lnN and lnDq in fully stocked stands. The maximum density was estimated using the 0.975th quantile [28] of the R function quantreg [29]. In addition, we also calculated the Reineke's Stand Density Index (SDI) for each plot (using our own empirical exponent 1.536, instead of Reineke's 1.605). The SDI has been used as a measure of density for predicting forest productivity [30]: where N is the number of individuals per area, Dq is the quadratic mean diameter (cm), and β is the allometric exponent that expresses the relation between tree size and number of trees estimated in Equation (1).

Forest Diversity
The variety of biodiversity indices is bewildering, and many are viewed with suspicion. Most diversity indices are entropies, not diversities, and their mathematical behavior often does not correspond to our intuitive concept of diversity [31][32][33]. However, there seems to be general agreement among statisticians that numbers equivalent [34], not entropy, should be the diversity measure of choice [35]. The first three Hill numbers depend on the weight that is assigned to rare species: q = 0 refers to species richness where rare species carry the same weight as common ones, q = 1 refers to the exponential of Shannon's entropy index, and q = 2 to the inverse of Simpson's concentration index. For example, if q = 1, the Shannon index H = pi × log(pi) is used, and the corresponding Hill effective number is N = e H' . Effective Hill numbers, instead of entropies, are used increasingly to characterize the taxonomic, phylogenetic, or functional diversity of a community. Obviously, estimates of Hill numbers are a function of the sampling effort.
Following Leinster and Cobbold [36], we may extend community diversity beyond the naive focus on a simple index. In this study, we used the number of tree species in plot (S) and effective Hill numbers to test the grouping of sites with different tree diversity. Obviously, the diversity of a forest ecosystem is not only characterized by species richness, but also by a set of attributes that portray the significant features of heterogeneity within a specific community [37]. The diversity profile of a forest community may be characterized by a number of attributes, including the distribution of species and sizes of individuals, and by the way in which these attributes are spatially mingled.

Modelling Forest Production
The assessment of forest production, i.e., the periodic annual volume increment per unit area and time (PAI), is fundamental to sustainable forest use [38]. PAI may decrease if the residual stand density diverges from the optimum [39], which may be caused by unsustainable harvesting [16,40]. Reduced production may also be expected in very dense forests where production may be offset by mortality and increasing competition [41]. To determine the relationship between forest site variables and PAI (defined here as the periodic annual volume increment of all live trees plus ingrowth (m 3 /ha/year)), we used a stepwise linear regression procedure. This method is useful for determining the relative importance of the independent variables and evaluating the order of importance of the variables [42]. The significance level of 5% was used for the selection of candidate predictive site variables in the model. These predictive variables were the measurements at the beginning of the selected five-year growth period. They included the number of trees per hectare (N), the Reineke's Stand-density Index (SDI); the mean plot height (Hmean, m), the Hill number as a diversity measure (Hill), the stand quadratic mean diameter (Dq, cm), the total basal area of the plot (G, m 2 ha −1 ), and the stand dominant height (H 0 , m). These predictive variables, along with their various combinations, were tested separately. Nonsignificant variables at 5% were removed in the modelling process. In addition, the relation of PAI and residual levels of G was studied to evaluate the effect of low forest densities on production in the study area.
Model evaluation is an important part of our analysis. We used the coefficient of determination (R 2 ) and the Residual Standard Error (RSE). In addition, the normality of the distribution of residuals and heteroscedasticity were evaluated using the Shapiro-Wilk test [43] and graphical analysis.

Stratification of Field Plots
The cluster analysis generated three groups of plots, which described distinctly different forest structures at Sierra Madre Occidental. Figure 2 presents four characteristics of the three groups of plots: (a) The relationship between number of trees and basal area; (b) the relationship between number of trees and dominant height; (c) the diameter distributions; and (d) the height distributions. For the same number of trees (the same "crowding"), forest densities differed greatly among groups. The scale and shape parameters of the Weibull function also differed considerably for diameter at breast height and total tree height.

Stratification of Field Plots
The cluster analysis generated three groups of plots, which described distinctly different forest structures at Sierra Madre Occidental. Figure 2 presents four characteristics of the three groups of plots: (a) The relationship between number of trees and basal area; (b) the relationship between number of trees and dominant height; (c) the diameter distributions; and (d) the height distributions. For the same number of trees (the same "crowding"), forest densities differed greatly among groups. The scale and shape parameters of the Weibull function also differed considerably for diameter at breast height and total tree height. The results of the Kruskal-Wallis Rank Sum Test, used to further discriminate among the three groups of field plots defined by the cluster process, indicate that the variables N, G, Dq, Ho, and PAI showed highly significant differences (p ≤ 0.01). Only the number of species (S) and Hill number did not show significant differences. This result show that at least one group was distinct from the others. For these cases, we applied Dunn's test to determine whether all three groups were different or whether only one group was different at the level of 5% (p ≤ 0.05; Figure 3). N, G, Dq, Ho, and PAI suggested significant differences for all of them, indicating that each group corresponded to a significantly different forest structure. The results of the Kruskal-Wallis Rank Sum Test, used to further discriminate among the three groups of field plots defined by the cluster process, indicate that the variables N, G, Dq, Ho, and PAI showed highly significant differences (p ≤ 0.01). Only the number of species (S) and Hill number did not show significant differences. This result show that at least one group was distinct from the others. For these cases, we applied Dunn's test to determine whether all three groups were different or whether only one group was different at the level of 5% (p ≤ 0.05; Figure 3). N, G, Dq, Ho, and PAI suggested significant differences for all of them, indicating that each group corresponded to a significantly different forest structure.  Table 1.

Maximum Forest Density
The relation between the number of individuals and the basal area, described by Reineke as a site density index that expresses the maximum tree occupancy in each area, was negative in all three density groups (Figure 4). The limiting density estimated for Durango´s forests was ln( ) = (118,434.6) − 1.536 × ln( ). The three groups differed from each other in stand density and squared diameter. The highest stand density was observed in Group 3, while the lowest was recorded in Group 1. The squared Figure 3. Forest conditions by class of density. Width represents relative density within each bin and the boxplot represents the quartiles. All the variables showed significant differences (p ≤ 0.05) according to Dunn's test, except for the number of species, which was not significantly different among groups. N, G, Dq, Ho, and PAI are defined in Table 1.

Maximum Forest Density
The relation between the number of individuals and the basal area, described by Reineke as a site density index that expresses the maximum tree occupancy in each area, was negative in all three density groups ( Figure 4). The limiting density estimated for Durango´s forests was ln(N max ) = ln(118, 434.6) − 1.536 × ln(Dq).  Table 1.

Maximum Forest Density
The relation between the number of individuals and the basal area, described by Reineke as a site density index that expresses the maximum tree occupancy in each area, was negative in all three density groups (Figure 4). The limiting density estimated for Durango´s forests was ln( ) = (118,434.6) − 1.536 × ln( ). The three groups differed from each other in stand density and squared diameter. The highest stand density was observed in Group 3, while the lowest was recorded in Group 1. The squared The three groups differed from each other in stand density and squared diameter. The highest stand density was observed in Group 3, while the lowest was recorded in Group 1. The squared diameter showed the opposite pattern. Dunn's test shows significant differences (P ≤ 0.01) among the three groups of plots (Figure 3).

Relationship between Forest Density, Diversity and Production
The stepwise linear regression identified the principal variables that determine the productivity in the multispecies and uneven-aged forests of the Sierra Madre Occidental. The equation selected to predict PAI is given below: where PAI is the periodic annual increment (m 3 /ha/year); N is the number of individuals per hectare; SDI is Reineke's Stand-density Index; Hmean is mean height; Hill is Hill number; and a, b, c, and d are the parameter estimates. Table 2 shows the parameter estimates, residual standard error (RSE), and coefficient of determination (R 2 ) of Equation (3). All parameters are significant (p ≤ 0.05) and the residuals were uniformly distributed. Interestingly, the coefficient d is negative, i.e., production is reduced considerably with increasing tree species diversity. RSE is the residual standard error and R 2 is the coefficient of determination.
According our PAI data, the residual density is a principal element affecting forest production. Figure 5 shows, as an example, the relation between basal area and production based on our data (PAI = 0.000788 × G1 3 − 0.0000130 × G1 4 , R 2 = 0.81). Any attempt to generalize or simplify this relation at high densities should be based on a careful analysis of additional variables, including species composition, structure, and climate. diameter showed the opposite pattern. Dunn's test shows significant differences (P ≤ 0.01) among the three groups of plots (Figure 3).

Relationship between Forest Density, Diversity and Production
The stepwise linear regression identified the principal variables that determine the productivity in the multispecies and uneven-aged forests of the Sierra Madre Occidental. The equation selected to predict PAI is given below: where is the periodic annual increment (m 3 /ha/year); N is the number of individuals per hectare; SDI is Reineke's Stand-density Index; Hmean is mean height; Hill is Hill number; and , , , and are the parameter estimates. Table 2 shows the parameter estimates, residual standard error (RSE), and coefficient of determination (R 2 ) of Equation (3). All parameters are significant (p ≤ 0.05) and the residuals were uniformly distributed. Interestingly, the coefficient d is negative, i.e., production is reduced considerably with increasing tree species diversity. RSE is the residual standard error and R 2 is the coefficient of determination.
According our PAI data, the residual density is a principal element affecting forest production. Figure 5 shows, as an example, the relation between basal area and production based on our data (PAI = 0.000788 × 1 − 0.0000130 × 1 , R 2 = 0.81). Any attempt to generalize or simplify this relation at high densities should be based on a careful analysis of additional variables, including species composition, structure, and climate. According to the relation between basal area and production, there is overwhelming evidence that production is significantly affected by forest density (Figure 5). This study shows that considerable loss in production can be expected if forest densities are reduced below 20 m 2 /ha. On the other hand, there are considerable potential gains in production if density is allowed to increase. The distribution of forest densities in our field plots and the density-production relation for Group 1 are shown in Figure 6. According to the relation between basal area and production, there is overwhelming evidence that production is significantly affected by forest density (Figure 5). This study shows that considerable loss in production can be expected if forest densities are reduced below 20 m 2 /ha. On the other hand, there are considerable potential gains in production if density is allowed to increase. The distribution of forest densities in our field plots and the density-production relation for Group 1 are shown in Figure 6. Assuming that the five-yearly basal area increment (∆G, m 2 /ha) may be estimated as a function of basal area (G) by the following equation (derived from this study, R 2 = 0.90, ∆G = 0.005287 × − 0.00006816 × ), it is possible to estimate the number of years required to raise an unproductive low density stand of 15 m 2 /ha to a more productive one of almost 25 m 2 /ha. The number of years required to reach and acceptable basal area of 25 m 2 /ha, assuming an unproductive initial basal area of 15 m 2 /ha, is shown in Table 3. At least 35 years without income from timber sales are needed to raise the unproductive basal area of 15 m 2 /ha to a more productive level of 25 m 2 /ha. Such examples demonstrate the dramatic effects of low forest densities for the communities which depend on the forest for income, and that low-density forest properties may need long recovery periods to attain acceptable levels of basal area.

Groups of Field Plots and Maximun Forest Density
The high natural forests in the Region Madrense in Durango State exhibited a range of complex structures and changing spatial patterns, and it was surprising that after more than 100 years of harvesting activities by the local communities, the tree species diversity showed still a high level of complexity. The results of this study did not reveal significant differences in tree diversity between well-defined field groups of plots that significantly differed in other site variables. However, the degree to which selective management may have caused changes in forest structure and species composition in the Region Madrense, which represents the most important ecotone of the SMO, has Assuming that the five-yearly basal area increment (∆G, m 2 /ha) may be estimated as a function of basal area (G) by the following equation (derived from this study, R 2 = 0.90, ∆G = 0.005287 × G 2 − 0.00006816 × G 3 ), it is possible to estimate the number of years required to raise an unproductive low density stand of 15 m 2 /ha to a more productive one of almost 25 m 2 /ha. The number of years required to reach and acceptable basal area of 25 m 2 /ha, assuming an unproductive initial basal area of 15 m 2 /ha, is shown in Table 3. At least 35 years without income from timber sales are needed to raise the unproductive basal area of 15 m 2 /ha to a more productive level of 25 m 2 /ha. Such examples demonstrate the dramatic effects of low forest densities for the communities which depend on the forest for income, and that low-density forest properties may need long recovery periods to attain acceptable levels of basal area.

Groups of Field Plots and Maximun Forest Density
The high natural forests in the Region Madrense in Durango State exhibited a range of complex structures and changing spatial patterns, and it was surprising that after more than 100 years of harvesting activities by the local communities, the tree species diversity showed still a high level of complexity. The results of this study did not reveal significant differences in tree diversity between well-defined field groups of plots that significantly differed in other site variables. However, the degree to which selective management may have caused changes in forest structure and species composition in the Region Madrense, which represents the most important ecotone of the SMO, has not been assessed in sufficient detail [17,24]. The preservation of species richness in these complex forest ecosystems requires new methods of management. Accordingly, Gadow et al. [3] presented new retention strategies for selectively managed natural forests. They found that methods that do not prescribe how much to harvest but specify the residual forest in terms of tree species and dimensions, the structure, density, and diversity remaining after the harvest, were especially relevant.
Several new permanent forest observational networks have been established since the turn of the 20 th century. Individual scientists have identified a real need for permanent monitoring that was not provided by the NFI's in their countries. An example of such an initiative is the Mexican permanent plot network, which provided essential information in this study. Another example is the Forest Observational Network of the Beijing Forestry University, which includes several very large field plots with mapped trees in the mixed deciduous forests of Northeastern China, in the northern pine-oak forests, and in the central and western Abies-Picea mountain forests of China [44]. Large permanent plots provide opportunities for evaluating effects of scale [28], while networks of regularly spaced small plots enable scientists to study environmental response gradients. Both types are useful for the evaluation of the effects of forest density and the degree of site occupancy on forest production.
In this study, three groups of plots, which described distinctly different forest structures at the SMO, were identified. These forest types showed significant differences in N, G, Dq, Ho, and PAI. This result may be an indication that differences in forest density are mainly the result of management. The differences might also be caused by forest gaps created by small-scale disturbances that occur in these forests, like windthrow, bark beetles, and forest fires. However, the site variables used in the Kruskal-Wallis Test probably had a high level of collinearity, and it is therefore necessary to use other statistical techniques like ordination methods (e.g., principal component analysis, Detrended Correspondence Analysis) to describe which of these structural variables are the underlying factors that cause plots dissimilarity. Several other studies found that site occupation is an important property in a multispecies natural forest, and density management is a key to the successful use of the unique forest resources of the Region Madrense [6,16,24]. On the other hand, dominant height has the advantage of being hardly influenced by stand management measures such as thinning [38], although it can be affected by very low or very high stand densities [45]. In addition, forest structural attributes, such as number of trees and average tree size, may affect biomass productivity in natural forests [46]. Production in multispecies natural forests normally depends on site conditions, structural attributes, species composition, and forest density, as well as the interaction between these factors [47][48][49].

Maximun Forest Density
The limiting density estimate for Durango in this study was ln(N max ) = ln(118, 434.6) − 1.536 × ln(Dq). Hernández et al. [50] obtained similar results for the natural forests in Hidalgo, Mexico, where the values of α and β were 105,550.7 and 1.535, respectively. Results of the maximum density analysis indicate that the relationship between the mean squared diameter and tree density well describes the phase of stand development and is related to self-thinning and crown closure dynamics of these species-rich forests [51].

Relationship between Forest Density and Production
Ambiguity may be caused by different interpretations of forest production: (i) Net production, i.e., the periodic annual volume (or biomass) increment of all live trees plus ingrowth and the volume of trees removed in thinnings; (ii) gross production, i.e., net production plus the volume of trees lost to mortality; and (iii) accretion, i.e., the volume increment of only those trees that survived during the growth period [52]. The differentiation between definitions of forest production is important, as each may result in different density-production patterns [52]. Confusions arising from previous work in forest density-production relationships are caused by the fact that the same density may include considerable differences in forest structure and diversity. The same forest density may include few large trees or many small ones. The same density may be found on poor or good sites. Accordingly, to be realistic, estimates of forest production in multi-age forests need to consider such different patterns, as exemplified by this study.
Different models for the relation between forest density and production have been used. Based on previous work, Allen and Burkhart [53] defined three theoretical relationships between forest production and density, as "increasing," "optimal," and "constant." They pointed out that this relation is usually uncertain for very high densities. The PAI model developed in the present study included the site variables N, SDI, Hmean, and Hill number in their formulation. However, the most important contribution of the variance, explained by Equation (3), came from the stand density index, which provides a measure of density that is independent of age and site quality [54]. Significant evidence has been presented that forest density influences production in plantation monocultures [55,56] and natural forests [16,24]. These results confirm earlier findings obtained with fewer observations about the potential theoretical relationships between forest density and production [38,57]. Our results show that considerable loss in production can be expected if basal areas of forests are reduced below 20 m 2 /ha. Schütz et al. [58] found that a residual basal area of 20-24 m 2 /ha was appropriate for successful beech regeneration in a forest characterized by a mosaic of irregular groups.

Conclusions
In this study, three groups of forest stands, which described distinctly different forest structures at the SMO, were identified. These forest types showed significant differences in the number of trees per hectare, basal area per hectare, quadratic mean diameter, dominant height, and periodic annual increment. Tree species diversity was similar between the three groups of stands. This result may be an indication that differences in forest density are mainly the result of forest management. According our study, the site variables N, SDI, Hmean, and Hill number were the the principal elements affecting forest production of these forests. One of the surprising results, confirming earlier studies (e.g., [24]), is the fact that the potential production in some areas was high, and exceeded 20 m 3 /ha/year provided a residual density of at least 30 m 2 /ha was maintained.
Although the results indicate that the continuous cover forestry management system used in Durango's community forests has not affected significantly tree diversity, studies in other regions of the world have shown that there is no unique and ideal diameter distribution in multispecies forests [59]. These forests may exhibit a great variety of structures with varying proportions of small-, medium-, and large-sized trees. Tree species compositions are often highly variable in natural forest communities. For this reason, it is necessary to develop new methods for ensuring sustainable production in natural forests based on forest density and particular species distribution.
Most of the local Mexican Ejidos and Comunidades in the study region have been able to withstand the pressure of commercial interests to convert the natural ecosystems to even-aged monocultures. An alternative method for controlling sustainable harvests in multispecies natural forests is an approach based on residual basal area. The general concept is based on the premise that a post-harvest residual density is defined and that the residual density is distributed over specific tree sizes and tree species [3]. Such a method would not only ensure that the productive potential can be sustained, but also that the main features of the natural system are maintained. Details of such a management system need to be worked out to meet the specific demands of the communal forests of the Sierra Madre Occidental.