Introduction

Growth is the result of the cumulative increase in dimensions (i.e., tree, stand, or plantation) over a given period. Growth accumulates and is expressed as yield (Telles-Antonio et al., 2022), and it is determined by the interaction between genetic characteristics and the biophysical factors of the site where individuals develop (Salas et al., 2016). Climate, soil, and topography influence phenotypic plasticity and the species-specific yield response, due to varying rates in physiological processes (Vogel, 2018).

In natural stands and forest plantations, the estimation of growth and yield projection is fundamental for proposing silvicultural practices that agree with tree conditions and site productivity (Burkhart & Tomé, 2012). Quantifying the effects of biotic and abiotic factors on growth and yield is essential for understanding structural variation and species composition over the rotation period, and for recommending appropriate silvicultural practices (Salas et al., 2016). In this regard, Brown (2007) in the United States found that slope and soil depth influence the dominant height increment in Pinus strobus L. plantations (<53 years), leading to the proposal of differentiated models for its estimation. In Spain, Cantero Amiano (2020) documented a variable yield response based on total height and diameter at breast height in stands of Pinus pinaster Aiton ssp. atlantica, because of the altitude and exposure conditions. Geographic orientation (aspect) promotes microclimatic variation that contributes to the creation of specific habitats for the development of certain species (Haire et al., 2022), as it affects soil biophysical processes and influences vegetation structure and species composition-factors that generally determine growth rate and yield (Gerding et al., 2006; Romero et al., 2014).

Growth and yield models are quantitative silvicultural tools that assist in describing the behavior of variables of interest at a given age or in designing silvicultural practices (Salas et al., 2016). Due to the lack of quantitative silvicultural information for the management of forest plantations in Michoacán, the objective of this study was to develop an explicit Growth and Yield System (TGYS) by incorporating site orientation (Aspect) as a random effect in plantations of Pinus montezumae Lamb. and Pinus pseudostrobus Lindl. in Nuevo San Juan Parangaricutiro (NSJP), Michoacán, Mexico. The main hypothesis is that growth conditions, influenced by topographic exposure and the environmental characteristics of the sampling sites, determine distinct growth dynamics and differentiated yield potential for each forest plantation.

Materials and Methods

The study area is located between coordinates 19° 17’ - 19° 30’ N and 102° 06’ - 102° 17’ W in an altitudinal range of 2 000 to 3 000 m. The topography is characterized by mountain ranges, with extrusive igneous geology and predominantly Andosol soils, characteristic of the Neovolcánica Tarasca physiographic subprovince. The climate is classified as Cw, with an average annual precipitation of 1 500 mm and mean temperature of 15 °C (Instituto Nacional de Estadística y Geografía [INEGI], 2017).

Data were collected from 14 P. montezumae plantations aged 7 to 34 years and 31 P. pseudostrobus plantations aged 7 to 39 years. Model fitting was based on the remeasurement (in 2021) of 154 temporary sampling plots of 400 m² (20 × 20 m), originally established as a chronosequence in 2019; of these, 58 plots corresponded to P. montezumae and 96 to P. pseudostrobus. In each sample unit, measurements were taken for diameter at breast height (DBH; cm) and total height (Ht; m) of each tree, which were then used to calculate average height (Ah; m), average dominant height of the five tallest trees (Hd; m), quadratic diameter (Dq; cm), basal area (BA; m²), stand density (N), and volume (V; m³).

Individual stem volume (V, m³) was estimated using the equations proposed by García-Espinoza et al. (2018, 2019) for P. montezumae [1] and P. pseudostrobus [2]. At each site, the geographic orientation of the terrain (Aspect) was recorded as zenithal (Z), northeast (NE), east (E), southeast (SE), southwest (SW), and northwest (NW).

Because TGYS consists of equations describing the relationships among key tree measurement variables and can simultaneously predict changes over time (Tamarit-Urias et al., 2014, 2019), the algebraic difference approach (ADA) was selected. This approach is based on a base model and an ADA equation that share parameter values, allowing for predictions or projections of the variables (Santiago-García et al., 2016).

A total of six compatible algebraic difference systems for Hd in their ADA form were obtained using Maplesoft^TM version 2015.1^® (Maplesoft, 2015), applying the method to define site quality curves with their respective site index (SI) classes at a base age or reference age (Eb) (Fierros-Mateo et al., 2017; Torres- Rojo & Magaña, 2001). Table 1 presents two growth hypotheses: (1) growth rate among sites is constant but the asymptote differs (anamorphic), and (2) growth rate varies among sites while the asymptote remains constant (polymorphic) (Hirigoyen et al., 2018).

Table 1. Compatible systems used for dominant height with the algebraic difference approach.

The system was selected based on the best statistical fitting values, where the estimation trend closely matches the observed data distribution, allowing for its readjustment into a projection equation for initial dominant height (Hd₁) as an additional common parameter. This approach assumes distinct growth potentials for each Aspect/Site combination (Wang et al., 2020).

Once the Hd system was selected, compatible systems for site density were fitted using the re-measured plot data. Basal area and volume were individually fitted to identify growth and yield differences by Aspect/ Site at plantation level, since performing the fitting at plot-level or per hectare reduces sensitivity due to the number of observations at each aggregation level. The compatible systems used and listed in Table 2 were obtained from Fierros-Mateo et al. (2017), González- Benecke et al. (2012), Hirigoyen et al. (2018), Parra- Piedra et al. (2017), Santiago-García et al. (2013, 2016), Tamarit-Urias et al. (2014, 2019), Torres-Rojo and Magaña (2001) and Uranga Valencia et al. (2018).

Table 2. Compatible systems to predict and project the current (age: A ₁ ) and future (age: A ₂ ) condition of the evaluated Pinus plantations.

The 19 compatible systems were statistically fitted using the R ^® 2022.07.0 Build > nlme > ML (Pinheiro, 2022) program (Pinheiro, 2022) through mixed-effects models (MEM). For some parameters, the combination of Aspect/Site was additively included as a groupinwg variable: Z/i, NE/i, E/i, SE/i, SW/i, and NW/i, where i is the sampling site. To obtain the system's random parameters, the marginal maximum likelihood approximation of the best linear unbiased predictor was used (Mehtätalo & Lappi, 2022).

The MEM analysis was chosen because it allows for the estimation of fixed parameters, which are global for the population, and random parameters, which enable specific contrast scenarios at each grouping level (Correa Morales & Salazar Uribe, 2016; Mehtätalo & Lappi, 2022). The fitting was performed by constructing a bivariate database (Pinheiro & Bates, 2000), as two measurements of the same variable were taken from the same individual at two different time periods; for example, the fitted structure of the compatible system [20] for V takes the form of system [22].

where: V ₁ and V ₂ = Vector of observations of initial and projected volume, or age zero and Ab in the i-th individual, respectively; f (BA _1i, Hd _1i, A _1i; p) = structure of the prediction model including the independent variables; g (V _1i, BA _1i, Hd _1i, A _1i; p) = structure of the projection model including the independent variables; ε _i ~ N(0, σ ²) = error term within the system, and, and p = parameter vector, in which the mixed-effects parameter (ω ₁) is specified, incorporating the grouping level additively (ω ₁ + ω _i), resulting in a value for the random parameter defined by group (ω _1i). The structure with the random effect defined as ω _1i ~ N(0, ɸ²):

To correct for heteroscedasticity and inconsistent error in the estimates, a power variance function was included within the MEM adjustment structure: v a r ε i j = v i 2 ∙ δ 1 ∙ 1 - f + f . Where var(ε _ij ) is the variance function evaluated using the variance covariate of the residuals from the predictor (v ₁), and δ ₁ refers to the variance function coefficient, which is specific for each level (δᵢ). The term f is an indicator variable with a value of 1 when used in the prediction equation and 0 when applied to the ADA expression (Mehtätalo & Lappi, 2022; Pinheiro & Bates, 2000).

Autocorrelation was corrected by incorporating an autoregressive moving average structure (ARMA(p,q)), into the system, where p and q represent the assigned autoregressive and moving average orders, respectively (Pinheiro & Bates, 2000; Quiñonez Barraza et al., 2018).

The best compatible system was selected by verifying that all parameters were significantly different from zero (p < 0.05) and by evaluating the likelihood ratio through ANOVA (using the nlme function) (Pinheiro & Bates, 2000). Subsequently, the systems were ranked based on the values of the coefficient of determination (R ² ), root mean square error, bias, log-likelihood (logLik), and the Akaike and Bayesian information criteria (Santiago-García et al., 2013).

Since the best statistical fits do not necessarily reflect the distribution trend of the observed data, the estimates from each system were compared against field data. Normality in the frequency of residuals was assumed to be approximately normal (given >2 000 observations), and homoscedasticity of the residuals was verified graphically (Correa Morales & Salazar Uribe, 2016).

Once the estimates were obtained using the fixed parameters for each component of the TGYS systems, the current annual increment (CAI) was determined by the difference between the projected value (y ₁) and the initial value (y ₂) of each tree measurement variable, divided by the difference between A ₂ and A ₁ ( C A I = y 2 - y 1 A 2 - A 1 ). The mean annual increment (MAI) was calculated by dividing each tree measurement variable (y _i) by age (MAI = y _i / A _i), while the technical rotation (tt) was defined as the point at which CAI = MAI (Torres-Rojo & Magaña, 2001).

Results

The model fitting analyses for each tree measurement variable, after selecting the parameter where the random effect was included and which yielded the best statistical indicators, are presented in Table 3. This table implicitly includes the likelihood ratio test values obtained through ANOVA.

Table 3. Parameters and goodness-of-fit statistics of the compatible systems for tree measurement variables in Pinus species.

The detailed breakdowns of each fitted model are provided in Appendix I, which may serve as reference seed values for future fitting of the proposed compatible systems for each variable. For both species (Appendix 1.1), the Schumacher expression [8] (anamorphic) was statistically the best for Hd; however, when comparing the estimates to the observed data in P. montezumae, the trend overestimated values at ages <15 years and again after 24 years. Therefore, expression [7] (polymorphic) was selected instead. The variability explained by the adjusted compatible systems for P. montezumae ranged from 76 to 92 %, with overall deviations of 1.82 m and a bias of -0.73 m. For P. pseudostrobus, the explained variability ranged from 85 to 98 %, with global deviations and bias values of 2.08 m and -1.08 m, respectively (Appendix 1.1). The best performance was obtained by including the random effect defined by the Aspect/Site conditions in the growth rate parameter (β ₁) of the model.

The polymorphic growth trend in Hd at Eb of 25 years for P. montezumae, with four site index labels of 10 m, 14 m, 18 m, and 22 m (Figure 1a) shows different maximum increments and technical rotations for each condition (Figure 1b). This indicates that management and silvicultural practices should be differentiated. For P. pseudostrobus, the growth trend is anamorphic (Figure 1c) at Eb of 20 years, with higher CAI between 14-15 years for SI labels of 10 m, 16 m, 22 m, and 28 m (Figure 1d).

The likelihood ratio test via ANOVA indicated no significant differences between the two density systems for either species (α < 0.5). However, compatible expression [9] yielded the highest R² and greatest parsimony, and, given the similarity in RMSE and bias values, it was selected as the best-performing model (Appendix 1.2).

This expression produced more conservative estimates for initial density (N ₁ ), with 1 306 and 1 221 trees·ha^-1 for P. montezumae and P. pseudostrobus, respectively. In both cases, the random effect is associated with the rate of density decline as the stand ages and trees increase in size (τ ₁). The annual mortality rate was estimated at 3.14 % and 3.35 % for each species, with stand densities at the base ages of 25 and 20 years of 607 and 640 trees∙ha^-1 for P. montezumae and P. pseudostrobus, respectively.

When fitting the systems for BA, it was found that parameter γ ₃ of system [13] was not significant in either species; this parameter is associated with the ratio of age to initial density (N₁). Additionally, for P. montezumae only, parameters γ ₃ and γ ₂-γ ₃ in the compatible expressions [12] and [16], respectively-both related to stand density-also lacked significance. Therefore, these expressions are not included in the statistical analysis presented in Appendix 1.3. For both species, system [15] was selected based on both statistical fit and estimation trends. Using this system, basal area estimates at the Eb of 25 and 20 years were 23.12 m²∙ha^-1 and 30.15 m²∙ha^-1 for P. montezumae and P. pseudostrobus, respectively, with CAI at age of 10 years of 1.87 m²∙ha^-1∙yr^-1 and 2.62 m²∙ha^-1∙yr^-1.

System [17] was identified as the most appropriate for estimating volume in both species, based on statistical fitting and ANOVA results (Appendix 1.4), as it includes a random effect in the parameter associated with A ₁ under each condition (ω₂). Using the compatible structures as part of a full TGYS system, projected volume yields at the established base ages were 109.47-223.12 m³∙ha^-1 for P. montezumae (25 years) and 150.72-284.86 m³∙ha^-1 for P. pseudostrobus (20 years) (Figure 2a and 2c, respectively). The estimates suggest that yields of up to 260.95 m³∙ha^-1 and 468.70 m³∙ha^-1 can be achieved at the maximum recorded ages of 35 years for P. montezumae and 40 years for P. pseudostrobus, respectively.

The maximum CAI in volume for P. montezumae ranged from 7.2 m³∙yr^-1 and 13.7 m³∙yr^-1, varying by SI due to the polymorphic nature of the dominant height growth curve. For site indices of 10 m, 14 m, 18 m, and 22 m, the age of maximum accumulation occurs at 11, 13, 16, and 19 years, respectively, with corresponding technical rotations of 20, 23, 26, and 31 years (Figure 2b). In contrast, for P. pseudostrobus, the anamorphic trend in volume indicates a maximum CAI at 15 years and a technical rotation of 22 years, with yields ranging from 12.0 m³∙yr^-1 to 23.1 m³∙yr^-1 for the lowest (SI = 10 m) and highest (SI = 28 m) site qualities, respectively (Figure 2d).

Annex I.5 shows the standard deviation and residual values of the parameter that included the random effects for each compatible system, the values of the parameters related to the ARMA(p,q) correlation structure, and the variance function (varPower) in the analysis using mixed-effects models. Moreover, the variance-covariance matrix values of the parameters were included to allow for the use of these systems through some type of calibration of the random parameters.

Discussion

The polymorphic and anamorphic dominant height growth trends satisfactorily described the observed variability for P. montezumae and P. pseudostrobus plantations, respectively. The polymorphic Hd trend for P. montezumae included in the TGYS differs from that reported by Zepeda Bautista and Acosta Mireles (2016) for natural stands of P. montezumae in San Juan Tetla, Puebla, Mexico, who proposed anamorphic growth trends for Hd. However, it agrees with the findings of Fierros-Mateo et al. (2017) for Pinus chiapensis (Martínez) Andresen plantations in Tlatlauquitepec, Puebla, Mexico. On the other hand, in Nuevo San Juan, Michoacán, Mexico, García-Espinoza et al. (2019) reported polymorphic Hd growth trends for P. pseudostrobus using stem analysis data, which contrasts with the trend found in the present study for the species based on remeasurement plots.

The differences in dominant height growth trends may be partially attributed to variations in the implementation of silvicultural practices, such as uneven pruning and brushing, the intensity of silvicultural treatments such as pre-thinning or thinning or the type of management such as total or clear-cutting. A similar pattern was observed in Pinus patula Schl. et Cham., where initial analyses indicated that anamorphic Hd curves provided the best fit for growth modeling (Santiago-García et al., 2013); however, in a subsequent update of the TGYS, a polymorphic trend yielded a better fit (Santiago-García et al., 2016).

The relationship between Hd and age, as expressed through SI, is the most important component of the TGYS (Santiago-García et al., 2016), as it enables a simplified and accurate classification of forest productivity and supports growth projections for the planning of silvicultural activities. This is because height is the variable least affected by stand density level (Tamarit-Urias et al., 2019; Torres-Rojo & Magaña, 2001).

The mortality rates for P. montezumae (3.14 %) and P. pseudostrobus (3.35 %) are similar to those reported by Santiago-García et al. (2013, 2016) for P. patula stands in managed forests (3.3-4.76 %). However, they differ from the findings of Fierros-Mateo et al. (2017), who reported a mortality rate of 1.7 % for P. chiapensis plantations, which may be attributed to species-specific traits, site conditions, or the initial planting density. When using these mortality equations, it is recommended to establish a minimum stand density per hectare at a given base age. Without this restriction, the projections could lead to a scenario of zero density, which is neither logical nor desirable (Fierros-Mateo et al., 2017; -García et al., 2016). For example, setting a harvest rotation based on investor objectives and market demand for specific log dimensions can help ensure more realistic outcomes.

One approach to expanding the systems in future model refinements could involve incorporating the initial planting density (N ₁ ), as proposed by Tamarit-Urias et al. (2019) for Tectona grandis L. f. in Campeche, Mexico. Alternatively, a fixed N ₁ value could be established using spacing-density combinations such as 2 500: 2 × 2 m; 1 667: 2 × 3 m; 1 111: 3 × 3; and 833: 4 × 3 m, similar to the approach used by Parra-Piedra et al. (2017), who set the initial parameter at 1 100 trees∙ha^-1 for P. patula in Zacualpan, Veracruz, Mexico.

The estimates of Hd, BA, and V generated by the species-specific TGYS structures exhibit a sigmoidal growth pattern, characterized by a defined inflection point at a particular age and a maximum horizontal asymptote-traits that are considered ideal for modeling forest growth variables (Kiviste et al., 2002). In addition, the proposed transition functions effectively capture and project changes in tree measurement variables of the TGYS from the current state of the stand (Santiago-García et al., 2013) and demonstrate analytical consistency among the functions employed (Zepeda Bautista & Acosta Mireles, 2016).

It should be noted that the technical rotations derived from the fixed parameters of the GYS for P. montezumae and P. pseudostrobus may vary depending on site quality, generally resulting in longer rotations for lower site indices (Fierros-Mateo et al., 2017). Typically, technical rotations are shortened when modified through silvicultural treatments such as thinning, which reduces the crossover age between CAI and MAI under lower stand densities or greater thinning intensity (González-Benecke et al., 2012; Santiago-García et al., 2016). Additionally, rotation age may vary due to growth trends influenced by species-specific responses to topographic orientation. In younger stands, silvicultural practices such as pruning, brushing, and thinning contribute to increased growth and yield, and therefore, the optimal rotation age for a plantation can be variable (Fierros-Mateo et al., 2017).

The trends in dominant height, mortality, growth, and yield observed in the study species and explained by Aspect/site conditions are consistent with the findings of Zepeda Bautista and Acosta Mireles (2016), who reported clear relationships between the productivity of Mexican conifers and the biophysical characteristics of the site. These results also agree with those of Gerding et al. (2006), who analyzed how soil, climate, and topographic features influence tree growth in temperate regions. Similarly, Romero et al. (2014) identified the optimal exposures for the development of Pinus cembroides Zucc. (south and southwest) and Pinus johannis M.-F. Robert-Passini (north and northeast) in San Luis Potosí, Mexico.

The proposed models of TGYS represent a reliable option for estimating growth and projecting yield, particularly when the values of the random parameters are calibrated for each specific condition, as suggested by Meng and Huang (2009) and Sirkiä et al. (2015) using different methods. Based on the results of this study, forest management can be adapted according to topographic orientation, since factors such as exposure create microclimatic variations that significantly affect growth and yield (Gerding et al., 2006). Moreover, assessing growth and yield by topographic orientation provides an opportunity to explain reductions in productivity or stand performance due to factors such as elevation, wind, adverse weather events, or the presence of pests, diseases, and forest fires (Cantero Amiano, 2020).

Finally, updating the TGYS is essential for the proper management of forest resources (Santiago-García et al., 2016), and therefore, these systems should be subject to periodic review. Moreover, integrating such models into a growth and yield simulation program can enhance the applicability of variable-specific structures, support the evaluation of diverse management scenarios for improved forest planning (Santiago-García et al., 2013), and assess the effects of silvicultural practices such as thinning (González-Benecke et al., 2012), particularly when considering site index or initial planting density (Tamarit-Urias et al., 2019). Likewise, the development of density management guides for these species would contribute to validating or complementing the results obtained from the proposed TGYS (Uranga Valencia et al., 2018).

Conclusions

The Timber Growth and Yield System (TGYS) proposed for Pinus montezumae and Pinus pseudostrobus plantations in Nuevo San Juan Parangaricutiro, Michoacán, Mexico, proved to be statistically valid and consistent with the biological growth trends of the species. Therefore, it can be reliably applied using the fixed parameters estimated for each component of the system. The TGYS allows for the projection of increment, growth, and yield scenarios based on topographic orientation and site conditions, which can support the planning of site-specific forest management strategies for P. montezumae and P. pseudostrobus. The estimated annual and average growth rates, and technical rotation ages for both Pinus species can also guide the planning of cultural practices such as weed control and pruning, as well as silvicultural treatments including thinning and final harvesting in these plantations.