Study area

The study area is located in the southwestern region of Durango state, Mexico, specifically in the forests of the Borbollones ejido, Pueblo Nuevo (Figure 1), which is located at geographical coordinates 23º 30’ and 24º 15’ N, and 105º 15’ and 105º 45’ W. The predominant vegetation type corresponds to mixed-species forests with mainly species of the Pinus and Quercus genus. The elevation ranges from 1 400 to 3 000 m. The main climate is temperate semi-cold with an annual precipitation regime ranging from 900 to 1 200 mm and the mean annual temperature ranges from 8 °C in the upper areas to 24 °C in the lower ones (^{García, 1981}).

Study data

The data come from 441 temporary circular sampling plots of 0.10 ha, collected in 2012 with a stratified random sampling design for timber management purposes, covering an area of 2 339.89 ha divided into 263 stands. With the information of all species recorded in each plot, the following stand variables were estimated: number of trees per hectare (N; trees·ha⁻¹), basal area (G; m²·ha⁻¹), quadratic mean diameter (d_g; cm), dominant height estimated as the average of the 100 trees with the largest diameter per hectare (H_d; m) (^{Assmann, 1970}) and total tree volume (V; m³·ha⁻¹). The species composition was 84 % for the genus Pinus (P. ayacahuite Ehrenb., P. chihuahuana Schiede ex Schltdl. & Cham., P. cooperi Blanco, P. douglasiana Martínez, P. durangensis Martínez., P. engelmannii Carr., P. herrerae Martínez., P. leiophylla Schlecht. & Cham., P. lumholtzii B. L. Rob. & Fernald., P. michoacana Martínez., P. oocarpa Shiede. and P. teocote Schlecht. & Cham.); 8 % for the genus Quercus (Q. candicans Neé, Q. crassifolia Ehrenb., Q. durifolia Seem., Q. eduardii Trel., Q. scytophylla Liebm., Q. obtusata Bonpl., Q. rugosa Neé, and Q. sideroxila Humb. & Bonpl.); and 8 % for species such as Juniperus deppeana Steud., Cupressus lusitanica Mill. and Pseudotsuga menziesii (Mirb.) Franco. Table 1 shows a summary of the main statistics of the variables involved in the construction of the DMD.

Equations used in the DMD

The basic structure of a DMD is composed of two allometric equations fitted from density data, and a stand thickness index that is usually used to characterize the silvicultural treatments to be performed. The first of the equations relates the diameter of the average tree with the density of the stand and a signal of its productivity. The second equation predicts the stand productivity volume based on the variables that define mean tree volume and density. In turn, the variables of the DMD axes are conditioned by the characterization index of the chosen thickness, which forces the density to be represented in terms of the number of trees per hectare. In this study we used the Hart-Becking index "HBI" (^{Wilson, 1946}), also known as the relative spacing index, which is defined as the ratio (expressed as a percentage) between the average distance between the stand trees and the dominant height. The HBI is usually used to define the upper and lower density limits, as well as the need, intensity and frequency of thinnings (^{Diéguez-Aranda et al., 2009}); in addition, it is independent of the age of the stand (^{Long, 1985}) and can be adapted for even- and uneven-aged forests (^{Gadow & Hui, 1999}), thus representing a good silvicultural tool for managing density in mixed-species stands (^{Barrio-Anta & Álvarez-González, 2005}).Starting from two coordinate axes corresponding to the dominant height (x axis) and the number of trees per hectare (y axis), the isolines of the HBI values and the dependent variables of the allometric relationships are represented: one that relates the quadratic mean diameter with the number of trees per hectare and the dominant height, and another one that relates the stand volume with the number of trees per hectare, the quadratic mean diameter and the dominant height as an indicator of the site quality. The expressions of the equations used in the construction of the DMD are the following (isolating density):

These relationships have silvicultural and biological consistency; that is, the value of the quadratic mean diameter will be conditioned by the silvicultural treatments (which determine the residuals of the number of trees per hectare) and by the forest productivity (characterized by the dominant height). On the other hand, stand volume is directly related to the basal area (whose value depends on the number of trees per hectare and the quadratic mean diameter) and dominant stand height.

For the HBI calculation, the mean distance between trees corresponding to a stand with a triangular tree distribution has been used, which is more appropriate to the current state of the studied stands (^{Castaño-Santamaría, Barrio-Anta, & Álvarez-Álvarez, 2013}). Considering the above equations, if N is isolated for each of them and different values of HBI, d_g and V are set, expressions can be obtained that allow estimation of the isolines or trajectories that are subsequently superimposed on a cartesian coordinate system with H_d on the x axis and N on the y axis, thus defining the DMD.

The allometric system defined by equations 2 and 3 has d_g as an instrumental variable; that is, d_g is the dependent variable of equation 2, while in equation 3 it is an independent variable. Therefore, both equations were fitted simultaneously to avoid error correlation. The parameters were fitted using NLS with the 'nlsLM' procedure of the minpack.lm package (^{Elzhov, Mullen, Spiess, & Bolker, 2013}) in R (^{R Core Team, 2014}).

The establishment of the self-thinning line or DMD upper limit is considered of significant importance, given its usefulness in simulating natural mortality or self-thinning and prescribing thinning; if this line is not considered, there is a high probability of mortality due to competition. Quantile regression offers the possibility of obtaining the maximum density line and creating different quantile lines for stands that exceed the values observed within the diagram. The maximum density line was fitted by non-linear quantile regression (^{Koenker & Bassett, 1978}) with a potential equation that relates the number of trees per hectare to the dominant height:

where, a ₁ and a ₂ are the parameters to be estimated by minimizing the following function:

The value of the quantile (τ) used in the fitting process was 95 %, i.e. the upper limit of self-thinning is only exceeded by 5 % of the stands. Parameters a ₁ and a ₂ were estimated using the “nlrq” function of the “quantreg” module (^{Koenker, 2015}) of program R (^{R Core Team, 2014}).

Evaluation of Allometric Relationships

The analysis of the fitting statistic capacity of the allometric relationships was based on the graphs of the residuals and on statistics such as coefficient of determination (R ²) and Root Mean Square Error (RMSE), whose mathematical expressions are as follows:

Allometric relationships

The estimates of the β _i parameters of equations 2 and 3 that define the DMD (Table 2) were different from zero at a significance level of 5 % (P < 0.05). The analysis of the residuals showed that the regression equations did not present anomalous tendencies that implied non-fulfillment of the starting hypotheses of normality, homogeneity of variance and independence of errors. By means of the simultaneous fitting of equations 2 and 3, the independent variables were able to explain close to 87 % of the observed variability of the dependent variables, a common situation for static stand models (^{López-Sánchez & Rodríguez-Soalleiro, 2009}), minimizing the value of RMSE to 2.3. These results are consistent with the ranges observed in the database and show a correct estimate of d_g and V with the form of the equations. The accuracy in the fitting of equations 1 and 2 is acceptable when comparing the R² values with those of other studies carried out on even-aged stands of one or two species (^{Barrio Anta & Álvarez González, 2005}; ^{Castedo-Dorado et al., 2009}; ^{Tewari & Álvarez-González, 2014}) and on stands with a mixture of species of the genera Pinus, Quercus, Cupressus, Alnus and Juniperus (^{Corral-Rivas et al., 2015}).

Maximum density limit

The estimated parameters of the self-thinning line (

) were obtained by quantile regression with an intercept of 52 304.1 and a slope of -1.4391 for the independent term and exponent of the dominant height (a₁ and a₂), respectively. Both parameters values were different from zero (Pr < 0.05) and the slope value of the maximum density line was very close to Yoda's "Law of -3/2” proposal (^{Yoda, Tatuo, Husato, & Kazuo, 1963}). Although this relationship may vary due to site factors (species, productivity, age and shade tolerance), in this study it was very consistent for species diversity and observed density range. In particular, this result supports the conclusion presented by ^{del Río, Montero, and Bravo (2001)} and ^{Comeau, White, Kerr, and Hale (2010)}, who state that the slope is not always close to the theoretical value and may differ significantly between species composition and site quality. Therefore, the slope of the self-thinning model should be estimated with the data for each species and region studied, since populations have a different mortality rate depending on density or growth habits (^{Santiago-García et al., 2013}). ^{Bi, Wan, and Turvey (2000)} and ^{Zhang, Bi, Gove, and Heath (2005)} indicate that the fitting of the maximum density line by NLS is sensitive to data selection and may produce a line with an inappropriate slope. By contrast, quantile regression solves the problem by minimizing an asymmetric function with the absolute error loss, focusing only on extreme data (^{Bi, Bruskin, & Smith, 2002}).

Density management diagram

In construction of the DMD, first the maximum density line was drawn using the equation

. Subsequently, the density isolines (equation 1) were superimposed, representing different growth zones on which the level of competition of a particular stand is evaluated. Mortality was presented approximately from the isoline corresponding to 28 % of the HBI; approaching 33 %, mortality was more evident, due to the effect of the competition of the trees for the site’s resources (Figure 2). On the other hand, in the DMD, the lower limit of the constant growth zone, for those stands that have a dominant height between 14 and 20 m, can be considered from the isoline corresponding to 40 % of the HBI. These results agree with the diagrams developed by ^{Corral-Rivas et al. (2015)} that define the lower limit of the free growth zone between 40 and 45 % and the maximum density (upper self-thinning limit) between 25 and 30 % of the HBI for the mixed-species forests of Durango. By contrast, ^{Barrio-Anta and Álvarez-González (2005)} found the self-thinning line for pure stands of Quercus robur L. from 14 % of the HBI and set the 36 % isoline as the lower limit. This finding suggests that the constant growth zone is between 33 and 37 % of the HBI; therefore, this interval can be considered the level where there is full occupation of the site which, consequently, is where the gross growth of the stand is maximized.

Figure 2 illustrates a case study for the use of the DMD in a management scheme, whose objective is to obtain mostly saw-timber at rotation age. In general, any density program in the diagram focuses on maximizing stand volume production; therefore, density management options can be represented graphically from the ordinate axis (assuming there is no natural mortality). Reducing density by thinning has no effect on the evolution of dominant height (abscissa axis), so the speed at which the mass moves (age of the dominant height) through the abscissa axis depends on the quality of the site or the growth (^{Reineke, 1933}). Total stand production can be obtained directly from any point on the diagram from the isolines representing the total volume. The volume extracted, which results from the various interventions in the form of thinnings, can be estimated as the difference between the volumes before and after the intervention. The sum of these volumes throughout the forest rotation can lead to the estimation of total production under a given silvicultural management regime.

Figure 2 shows a scheme of silvicultural treatments (routes between points A-B, C-D and E-F) in a particular stand with the purpose of obtaining saw-timber at the final harvest (point G of the diagram), which is defined by a dominant height close to 20 m and a quadratic mean diameter of 20 cm, which will give an average volume of 180 m³·ha⁻¹. The scheme consists of three interventions (the first corresponds to a commercial thinning) arranged over time by means of the biological criterion of growth in dominant height of more than five years (routes between points B-C, D-E, F-G), a reasonable value for species that grow in good site conditions; to reach the final objective (point G), the scheme is presented in a stairway-shaped line (dotted line). The vertical segments represent the thinning cuts, while the horizontal ones reflect the intervals between thinnings in a 10-year cutting cycle, assuming that after thinning there will be no mortality due to competition among individuals; however, it should be considered that there is always mortality due to pests, diseases or meteorological factors, which is not considered in the DMD. Since the time variable to define the interval between cuttings is not shown in the diagram, the dominant height for the cutting cycle was estimated using the equation developed by ^{Corral-Rivas, Álvarez-González, Ruíz-González, and Gadow (2004)}; in addition, the number of silvicultural treatments applied to the even-aged stands for a 60-year forest rotation period was considered. According to the diagram, the first intervention (point B) should be carried out when the stand reaches a quadratic mean diameter of 11 cm and 1 200 trees·ha⁻¹ with an HBI of 25 % which, on average, represents a volume of 117 m³·ha⁻¹. The intermediate cutting segments were represented parallel to the ordinate axis under the assumption that the thinnings at the bottom do not affect stand dominant height (^{López-Sánchez & Rodríguez-Soalleiro, 2009}). As for production, this can be obtained at any point of the diagram using the stand volume isolines. It is also possible to estimate the wood volume, according to the industrial destination, from the commercial volume equations fitted for the main species of the study area (^{Pérez-Cruzado, Merino, & Rodríguez-Soalleiro, 2011}). The interpretation of the thinning program can be seen in greater detail in Table 3.