Modeling tree diversity, stand structure and productivity of northern temperate coniferous forests of Mexico

Study area

This research was conducted in the main core of the Central Portion of the Western Sierra Madre Mountain Range of Durango, Mexico. The region is characterized by a cold-temperate climate, with annual average rainfall and temperature of approximately 1,000 mm and 14 °C, respectively. Soils are classified as Litosols, Rendzins, Leptosols, and Regosols. Soils are shallow (<30 cm in depth), rich in acidic bases, and high in organic matter content. Unevenly aged, mixed coniferous forests with Pinus cooperi, Pinus durangensis, Pinus leiophylla, Pinus teocote, Quercus sideroxylla, Arbutus xalepensis, Arbutus chiapensis, Juniperus spp, Alnus accuminata, and other less frequent tree species distribute along this mountain range. Permanent sampling plots are located at an altitude of 2,450 m above sea level.

Permanent plots and silvicultural treatment

In this forest range, 36 permanent sampling plots were demarcated in 1967 in the municipality of San Dimas, Durango, Mexico. Each permanent plot was 100 m × 100 m quadrat with a measurement sub-quadrat of 50 m × 50 m located at the center for continuous inventory of dasometric parameters. Plot selection was conducted following a complete randomized design with six treatments and six replicates. Silvicultural treatments consisted on removing basal area at different intensities; (0) 0% or control, (1) 20%, (2) 40%, (3) 60%, (4) 80%, and (5) 100% or clear-cutting. Dasometric data recorded during the forest inventories of 1982, 1993, and 2004 were used to derive parameters of stand productivity, tree diversity, and stand structural complexity.

Forest productivity

Forest productivity was estimated for below ground biomass (coarse roots, BGB) and aboveground (stems, branches, and leaves) standing biomass (AGB), hereon denoted productivity, for three periods (1982–1993, 1993–2004, and 1982–2004) using allometric equations developed by Návar (2009) for large pine trees (AGBp = 0.0726D^2.4459) and for oak trees (AGBq = 0.0768D^2.4416). Allometric equations estimated standing AGB for all coniferous and all broadleaf tree species. Assessments were compared to independent calculations of biomass by the classical physics equation M = VP; where M = biomass; V = volume, and P = wood specific gravity, wsg. The Schumacher volume equation derived previously by Contreras-Aviña & Návar-Cháidez (2002) was used to evaluate timber volume. According to a literature review, the average wsg for pines was 0.43, and for oaks was 0.60 g cm⁻³ (Davalos, Wangaard & Echenique-Manriquez, 1977; Compean-Guzman, 1996). This comparison tested the consistency of both AGB estimates . The Clark et al. (2001) technique was used to estimate stand productivity following method 2. Using the allometric equations, total AGB was estimated for 1982, 1993, and 2004. Forest productivity was derived from (AGB1993-AGB1982)/11 for the first period; (AGB2004-AGB1993)/11 for the second; and (AGB2004-AGB1982)/22 for the full period of observations. The same procedure was employed to evaluate coarse root biomass productivity, but using the coarse root equation (BGB = 0.016D^2.668) reported by Návar (2009).

Tree diversity

Tree diversity was parameterized by fitting nine diversity indices for data recorded during each of the three forest inventories for each quadrat: three diversity indices based on species richness, one based on abundance, and four based on the proportional abundance of species. The diversity indices employed were species richness (S), Margaleff (Mg), Menhinick (Mn), Shannon Weiner (SW), Brillouin (Br), Simpson (Si), McIntosh (Mc), and Berger–Parker (B-P).

Four models fitted the abundance-diversity relationship for each plot for each of the three forest inventories to explain partition of limiting resources within the species community, allowing the drawing of conclusions regarding the successional processes that shape the productivity-diversity relationship. The equations describing the geometric series, the log series, the truncated lognormal distribution, and the broken stick model reported in McArthur (1957) and Magurran (2004) were fitted for each quadrat for each of the three forest inventories.

Stand structure

Various measures of stand structure have been proposed for guiding forest management (Staudhammer & Lemay, 2001). In this research, the three-parameter Weibull distribution for the variables top height, H, and diameter at breast height, D, depicted well the tree-size distribution that portrays the three-dimensional structure of forests (Haan, 2003; Hahn & Shapiro, 1967). The conventional procedure of moments in software reported by Návar & Contreras (2000) evaluated the shape, α, scale, β, and location, ε, parameters of the Weibull density function for D and H for each forest plot for each of the three forest inventories.

Stand attributes, diversity, and productivity

Linear and quadratic second-degree polynomial regressions fitted AGB productivity, P, to predict diversity, S, [(e.g.,  S = a + b(P) + c(P)²] for each of the two biomass components as well as for each of the three time periods using data at the quadrat scale. Comparisons of the coefficients of determination, r², and the significance of the model determined the type of regression that best describes the tendency of these relationships. These statistics test the appropriateness of the Unimodal, U-shaped and monotonic models (Mittelbach et al., 2001). In addition, I tested whether the quadratic relationship reached a maximum or minimum within the observed range of forest productivity values. This test examines the null hypothesis that c = 0, and checks whether the data supports an intermediate maximum or a minimum number of species within the range of observed productivity. This test also provides information on whether the quadratic is a better fit than the monotonic increasing or decreasing diversity with productivity (Leibold, 1999). Quadratic relationships that showed a maximum, negative c coefficient were classified as hump-shaped, while those that showed a minimum, negative b and positive c coefficients, were classified as U-shaped. Relationships without a statistical significant maximum (hump-shaped) or a minimum (U-shaped) were classified as linear models since transformations can linearize any monotonic relationship. All nine diversity indices were regressed against below and aboveground standing productivity. Therefore, a total of 54 regressions were tested for linearity or non-linearity (9 diversity indices and 2 productivity component) for all 36 quadrats for each of the three time periods between forest inventories.

The growth and yield model

Pines (P. cooperi) and oaks (Q. sideroxylla) are dominant tree species in these forest plots. Two types of timber growth and yield models were used to test the binary (pines and oaks) effect of tree diversity on forest productivity (Clutter et al., 1983; Vanclay, 1994). The first type of growth and yield model was constructed with the typical variables site index, basal area, the average age of pine trees, and the diversity component given by the density of oak trees/total stand density. The effect of the diversity component on timber growth and yield was tested by establishing the null hypothesis that the coefficient = 0. The second growth and yield model was developed individually for the six stands treated with clear-cutting (stands with no oak trees) and the rest of the forest stands (stands with oak trees of differential density). Timber volume projections over time by both types of models tested the null hypothesis that timber growth and yield are equal over time between mixed stands and mono-specific stands with only pine trees.

Stand structure complexity and productivity

The null hypothesis of the effect of canopy structure on forest productivity was tested by fitting regression equations between the Weibull parameters (α, β, and ε) of the H and D stand parameters and forest productivity for both time periods. The probability of the regression coefficients tested e.g., Ho:B1 = 0 or Ha:B1 ≠ 0. If Ho is rejected, then structural features represented by the Weibull parameters of the canopy explain part of the variation in forest productivity. In addition, statistical relationships between H and D and forest productivity were developed.

Forest stand attributes

Average stand timber basal area, volume, density, and root and aboveground biomass productivity components increased over time, unlike D and H, which oscillated erratically between time periods. Tree mortality, in-growth or recruitment, and likely errors in measurements of D and H accounted for part of the erratic averages of D and H. Pine trees dominated forest canopy in all 36 measured forest stands as basal area, stand density, timber volume, root biomass and aboveground standing biomass recorded larger statistics in pine, as compared to oak trees (Table 1).

Allometric and physical equations provided compatible AGB biomass assessments for oak and pine trees. Therefore, the conventional procedure of biomass estimation using allometric equations for individual trees was used to evaluate stand productivity. Total AGB productivity had an average (± confidence interval) of 3.09 ± 0.72 Mg ha⁻¹, 2.55 ± 0.77 Mg ha⁻¹, and 2.82 ± 0.60 Mg ha⁻¹ for the periods of 1982–1993, 1993–2004, and 1982–2004, respectively. Coarse root biomass productivity had an average (± confidence interval) of 0.17 ± 0.001 Mg ha⁻¹; 0.01 ± 0.0001 Mg ha⁻¹; and 0.08 ± 0.001 Mg ha⁻¹, for the time intervals of 1993–1982, 2004–1993, and 2004–1982, respectively. The drought spell of the 1990s (1989–2001) recorded and reported by Návar (2015) may have reduced total AGB productivity from 3.09 to 2.55 Mg ha⁻¹ y⁻¹, as well as coarse root biomass productivity from 0.17 to 0.08 Mg ha⁻¹y⁻¹, between the time periods of 1993–1982 and 2004–1993, respectively.

Removal of the basal area played a significant role by enhancing below and aboveground productivity (p ≤ 0.05; Table 2). There was an increasing tendency towards below and aboveground biomass productivity with basal area removal, with coefficients of determination (r²) of 0.36, 0.04, and 0.23; and 0.48, 0.25, and 0.44 for 1982–1993, 1993–2004, and 1982–2004, respectively. According to the tendency of these relationships, BGB and AGB productivity was least in control treatments, accelerated in a non-linear fashion with basal area removal, reaching the highest values in stands treated with 100% basal area removal. Growth of juvenile forests and most importantly the recruitment of young trees in forest openings partially explained the augment of productivity with increased basal area removal.

Tree diversity

P. cooperi regenerated in stands with 100% of the basal area removed, as well as several other stands with smaller openings. Q. sideroxylla also colonized canopy openings caused by tree mortality. Assemblages of seral tree species recorded during inventories were composed by P. teocote, P. leiophylla, P. ayacahuite, Arbutus spp, Alnus spp, Juniperus spp, Pseudotsuga spp, and Q. sideroxylla. Thus, oaks colonized forest openings and regenerated well under the canopy of pine trees in treated stands, as well. Tree density and species richness increased over time. Recruitment of individuals of secondary pine and oak species (P. ayacahuite, P. leiophylla, and Q. sideroxylla) enriched stand density and tree diversity. Most secondary species of succession are shade-tolerant and colonize well-stocked forest stands, as observed by mortality and canopy overlap between trees of the dominant pine species, P. cooperi. In control plots, Pseudotsuga menziesii colonized overstocked stands from 1993 to 2004. Species richness was accelerating during this period although the rate had decelerated since species richness, S, was increasing on the average from eight to 11 species per plot. Stand density was still increasing although the rate has decelerated notoriously causing the diversity indices to oscillate erratically over time (Table 3).

Stand density and species richness increased over time in all treated quadrats, except for density in stands with 100% basal area removal (Table 4). This tendency resulted from tree regeneration in gaps originated by basal area removal and the recruitment of shade-tolerant species under the canopy of remnant trees. In plots treated with 100% basal area removal, after reaching a maximum density of 1800 trees ha⁻¹ with an average D of 18. one cm in 1982, intrinsic competition between individuals of P. cooperi caused considerably tree mortality, reducing density to 1300 trees ha⁻¹ with an average D of 18. five cm in 1993. Shade-tolerant tree species P. ayacahuite, Q. sideroxylla, P. teocote, Junniperus spp, Alnus spp, Arbutus spp, and P. leiophylla colonized overstocked stands accelerating stand density to 1400 trees ha⁻¹ during 2004 through the reduction of P. cooperi trees. Intra-specific competition of P. cooperi reduced stand density while recruitment of shade-tolerant species increased stand density, as well as tree diversity.

At the plot scale, within a time interval, the diversity indices of species richness, Shannon & Weiner, Brillouin, Simpson, McIntosh, and Berger-Parker had a tendency to diminish over time. This was the result of increasing stand density and the decline of the rate of species richness with increasing basal area removal. In general, removing 100% of the basal area resulted in the smallest diversity indices (Table 4).

Most diversity indices however increased over time within each plot for all periods of recording data (Table 4). This tendency is noted primarily for all diversity indices in stands treated with the removal of 100% of the basal area. An increasing number of seral species colonized stands and the strong intra-specific competition of individuals that reduced the density of the pioneer species P. cooperi explained this tendency. Overcrowding processes of ecological competition reduce the density of dominant and increase the density of secondary species of succession. Although plots treated with total clearing increased species richness and revived stand density due to the regeneration of seral species, the diversity indices have not yet recovered when contrasted with those of other plots.

Diversity-Abundance models

The effect of basal area removal on the diversity—abundance structure is demonstrated by the number of significant null hypotheses accepted (p ≥ 0.05) for the geometric, log series, the lognormal, and broken stick models. Although there is no overall tendency for all of the accepted null hypotheses, it is clear that plots treated with 100% basal area removal have not yet recovered the abundance-diversity structure, since these stands had only 25% of all potential accepted null hypotheses. The rest of the stands had at least 50% of the accepted null hypotheses of the four fitted diversity-abundance models.

Diversity-abundance of most quadrats is developing quite effectively since the geometric series (a model adapted well to pioneer stages of succession) has the least number of null hypotheses accepted. Fitting tests accepted at least 50% of the null hypotheses of the broken stick model and 25% of the log series and lognormal model; three diversity-abundance models describing well the advanced stages of succession. That is, the differential abundance of most tree species is being reduced over time in these forests.

The canopy structure

The Weibull distribution fitted well D and H distributions for 92.5% of plots (probability of Kolmogorov-Smirnoff ≥ 0.05). Parameters describing the canopy structure showed the H structure developed over time from a nearly J inverse (a = 1.89 for 1982 and a = 2.19 for 2004) to a bell-shaped distribution with an increasing dispersion over time (b = 11.63 for 1982 and b = 14.44 for 2004). The smallest top height of trees was in the range of 2 to 3 m. Although Weibull parameters varied over time, the mean H (standard deviation) remained quite constant throughout the measurements. It was 9.01 m (2.61 m) for 1982, 8.88 m (2.56 m) for 1993, and 8.91 m (2.62 m) for 2004. Using the Weibull probability density function, there was a 74% and only a 13% probability of finding trees with 5 ≤ H ≤ 20 and H ≥ 20 m, respectively. Quercus spp, Arbutus spp, and Alnus spp occupied the lowest strata, H ≤ 10 m, followed by oaks in the middle strata with an average top height of 11.3 m (±4.3 m), and pine trees in the upper strata with an average top height of 13.0 m (±5.7 m). In forest openings, pioneer pine trees (P. cooperi) colonized well these places. Shade-tolerant pine trees (P. ayacahuite and P. leiophylla) established well under the canopy of other pioneer pine species.

Diameter structures also developed towards a bell-shaped, symmetric distribution but the rate of change was quite small (note the shape parameters of the Weibull distribution). Diameter structure contracted over time, in contrast to top height structures. Using the Weibull density function, there was a 40% and only a 28.7% probability of finding trees with 10 ≤ D ≤ 20 cm and D >20 cm, respectively. The probability of finding trees with D >30, 40, and 50 cm reduces to 10, 2, and 1%, respectively. Oaks recorded an average diameter and basal area of 17.0 cm (11. five cm) and 11.9 ± 7.7 m² ha⁻¹ while pine trees, characterized by P. cooperi, had an average diameter of 17.7 ± 10. 1 cm and basal area of 16.1 ± 7.7 m² ha⁻¹.

Trees in plots treated with 100% basal area removal reached a bell-shaped distribution for H and D in the last inventory. The variance of H and D was the smallest in this treatment, mimicking a well-stocked forest plantation. In contrast, untreated forest stands developed towards a J-inverse diameter distribution and to a bell-shaped H distribution. Mortality of large trees, recruitment of pioneer tree species in openings, recruitment of shade-tolerant tree species below canopies of large trees, and other unaccounted factors such as random tree mortality by competition enhanced canopy structure complexity. This tendency was also consistent for stands with 20% and 40% basal area removal.

The relationship between tree diversity and forest productivity

Most relationships between diversity indices and productivity did not present statistical significance (p > 0.05), with the exception of the monotonic increment of the Shannon Weiner and Brillouin indices with increasing productivity (p ≤ 0.0024). Correlations were positive but following a logarithmic tendency, with small coefficients of determination, r² < 0.10. All diversity indices increased, as productivity did, but the variation was so large that resulted in a lack of statistical significance. This was probably due to the fact that between 7 and 12 tree species colonized these forests with only two species dominating stand density (Q. sideroxylla and P. cooperi). Only a few individuals of other less common tree species (Alnus spp, Arbutus spp, and Pseudotsuga spp) were inventoried in few quadrats.

Timber growth and yield model

The relative density of oaks increased significantly timber growth (p ≤ 0.05) highlighting the relative importance of tree diversity on stand productivity (Table 5). Timber volume (m³ ha⁻¹) was a function of basal area (BA), site index (SI), the average age of pine trees (t), and oak density/stand density (IQ). However, IQ was better correlated to timber growth than the conventional SI derived from the Schumacher model fitted to the relationship of the top height of dominant pine trees and age. A simple representation of the two timber growth models, one for only pine forests and the second for mixed pine-oak stands is depicted in Fig. 1. Timber growth reached higher yields in mixed stands in contrast to mono-specific pine stands, although the effect appeared to be temporally limited. For the range of observations, timber growth was larger in mixed stands with the maximum deviance of 58 m³ ha⁻¹, reached at 25 years of stand projection. After 65 years of stand development, timber growth appeared to be higher in mono-specific pine stands.

Canopy structure and stand productivity

In highly productive plots, the height of most trees was restricted 5 ≤ H ≤ 15 m, and the horizontal distribution of trees appeared to depict well a J–inverse shaped distribution. In the least productive plots, on the contrary, the H distribution increased its variation (2 ≤ H ≤ 25 m) and was less balanced; the diameter distribution returned to a right-skewed distribution, indicating that the H distribution of trees appeared to be more important in determining AGB productivity.

Conclusions and recommendations

Observations derived from this research reveal forest productivity can be optimized with silvicultural prescriptions that must aim to balance the mixture of pines and oaks, as well as to conserve the structural complexity in the vertical and horizontal dimensions of trees. Tree diversity should follow a well-balanced combination of species richness and abundance that resembles intermediate rather than a pioneer stage of tree succession, mimicking probably diversity-abundance structures predicted by broken stick models. Top height structural complexity should be skewed with high variation (between 2.5 to 30 m), while diameter distributions skewed to the right. These recommendations are supported by the statistical relationships of forest productivity, tree diversity, and canopy structural complexity. These silvicultural prescriptions could preferentially buffer against contemporaneous perturbations associated with coupled heat waves and dry spells.