The study was conducted at the Intensive Carbon Monitoring Site Atopixco (SMIC-Atopixco) located in Zacualtipán de Ángeles, Hidalgo, Mexico. The SMIC-Atopixco is part of the Mexican Network of Intensive Carbon Monitoring Sites (RED Mex-SMIC) (Ángeles-Pérez et al., 2015) comprising a 900 ha (9 km²) polygon divided into nine 100 ha squares. It is located between the extreme coordinates 20°37´49.78” and 20°35’18.74” N and 98°37’51.01” and 98°34´2.71” W (Pérez-Vázquez et al., 2021), and mostly consists of portions of four ejidos of the municipality of Zacualtipán de Ángeles: La Mojonera, Atopixco, El Reparo and Tzincoatán (Ortiz-Reyes et al., 2015).

The climate is predominantly temperate humid (C(m)) followed by temperate sub-humid (C(w2)), with a mean annual temperature between 12 and 18 °C and mean annual precipitation between 700 and 2500 mm (Cruz-Leyva et al., 2010). In the managed forests, there are many species such as Pinus patula Schiede ex Schltdl. & Cham. and Pinus teocote Schiede ex Schltdl. & Cham., and in areas of natural vegetation, there are species such as Quercus excelsa Liebm., Quercus obtusata Bonpl., Quercus crassipes Humb. & Bonpl., Q. rugosa Née, Alnus sp., Clethra alcoceri Moc. & Sessé ex DC., Crataegus pubescens C. Presl., and Arbutus xalapensis Kunth (Cruz-Ruiz, 2004).

A flight was conducted in May 2013 to record LiDAR data in the study area. The Rielg VQ-480 sensor was used at an altitude of 397 m, pulse frequency of 200 kHz, scan angle of ±15° and flight line overlap of 50% (Ortiz-Reyes et al., 2015; Galeote-Leyva et al., 2022). The point cloud had an average density of 23 points m². LiDAR data processing was carried out using FUSION/LDV^® software version 4.50 (McGaughey, 2023) and FugroViewer version 3.5 (Fugro Geospatial Services, 2023) was used for three-dimensional visualization.

The LiDAR points positioned on the terrain were used to generate the digital elevation model (DEM) with the GridSurfaceCreate command of FUSION at a resolution of 0.8 m. Subsequently, using the ClipData command, the heights of the recorded returns were normalized, a process that consists of subtracting the DEM from the point cloud. A new point cloud was generated and clipped for the permanent sampling subunits (400 m² circular sites) present at SMIC-Atopixco using the PolyClipData command, generating 160 new files (one for each site). Then, for each file, LiDAR elevation metrics were obtained with the CloudMetrics command, which were related to the DH recorded in the field to identify the relevant predictor variables by means of a Pearson correlation matrix in R software (R development Core Team, 2023). Fifteen height percentile predictor variables (1st to 99th height percentile) and two descriptive variables (minimum and maximum height) were used. The stepwise regression methodology was tested as a complement to define the variable that best estimated DH in the linear model (L1), with the 99th height percentile (p99) being the variable that provided the best fit. The p99 was integrated into non-linear models (exponential (L2) and potential(L3)) to predict DH at site level (Table 1).

For the total area of the SMIC-Atopixco, LiDAR metrics were obtained in a grid of 20 m equidistant points using the GridMetrics command. From these, and to be consistent with the size of the sampling subunits, a vector mesh of square polygons with a surface area of 400 m² was generated in QGis version 3.28.11 (QGIS Development Team, 2023) and adjusted so that each point was in the center of each polygon. The LiDAR metric information associated with the points was copied to each polygon using the join attributes by location tool.

The SMIC-Atopixco is composed of 40 permanent sampling units, which are systematically distributed (Figure 1). Each of these units contains four 400 m² sampling subunits (sites) arranged in an inverted “Y” shape. This design is similar to the one used in the National Forest and Soil Inventory (INFyS) (CONAFOR, 2012) and follows the sampling scheme proposed by Hollinger (2008). Site data were processed assuming spatial independence among the study variables.

The data used in this study come from four inventory re-evaluations carried out from the initial measurement in 2013 (2014, 2015, 2019 and 2022) in the 160 permanent sampling sites of the SMIC-Atopixco. The dendrometric variables recorded at tree level include total height (H, m) and diameter at breast height (DBH, cm), while age (in years) was obtained at stand level according to the forest management program.

An audit was carried out on the database obtained from the forest inventory in the field, where the species P. patula was filtered and six individuals with outliers in the DBH-H relationship were identified, which were discarded from the analysis because they represented only 0.22%. The 2013 measurement, corresponding to the LiDAR scanning date, was used in a dynamic table in Microsoft Excel^® to obtain the DH of the 400m² sites (the average of the five tallest trees according to the suggestion of Bengoa (1999)), which were incorporated into the database with the LiDAR metrics per site.

From the total observations, a sample size of 329 (12.5%) trees was selected at 95% reliability to sort their information under the non-overlapping pairs method and with this, nine DH growth models were fitted under the ADA approach (Clutter et al., 1983; Cieszewski & Bailey, 2000) (Table 2). Each of the models was structured in the form DH₂, = f E 1 , E 2 , DH 1 β i , where DH 1 is the DH to be predicted at age E 2 ; β i are the parameters to be estimated DH 1 and E 1 and are the initial conditions of height and age, respectively.

Nine SI models were generated from the Hossfeld IV, Chapman-Richards and Korf base models using the ADA approach (Table 2). The base equations were rewritten in two different state conditions. To generate families of anamorphic SI curves, the parameter determining the asymptote in the base equation was cleared and substituted into the equation generated in the second state, whereby the asymptote is considered to be implicit, whereas, to generate families of polymorphic SI curves, the remaining parameters were isolated and substituted into the equation generated in the second state (Clutter et al., 1983). The nine DH growth models were fitted to the base age of 40 years as it was the rotation age assigned in the forest management program (FMP) for P. patula.

The SI models were fitted in R software (R development Core Team, 2023) with the nonlinear least squares (nls) procedure. The goodness-of-fit statistics analyzed were based on results of quantitative analyses; the highest value of the adjusted coefficient of determination R adj 2 , the lowest value of the root mean square error (RMSE), the lowest value of the Akaike information criterion (AIC), the significance of the parameters, and the distribution of the residuals were considered. Since the quality of the fit does not necessarily reflect the precision of the estimates of the observed data (Kozak & Kozak, 2003), a graphical comparison of the SI curves was carried out with respect to the observed data and the model that best reflected reality was selected. Five SI labels corresponding to 16, 21, 26, 31 and 36 m were used and classified as: high SI (>26m), medium SI (26 m) and low SI (<21 m).

The selected model was adjusted using the nonlinear mixed effects model (MEM) (with the grouping factor by plots) technique according to the suggested by Fang & Bailey (2001), for this, the following model was created: H D 2 = f A i β + β i b i , t i j + ε i j , where HD₂ represents the estimated dominant height in the i-th age of measurement; f() is the equation in ADA to be selected in the previous phase; A_i as the design matrix of size r×p for the fixed effects parameters; β is a vector of size p×1 containing the fixed parameters; B_i is the design matrix of size r×q for the specific parameter (random effect); b_i is a vector of size q×1 containing the random effects associated with the i-th plot. t_ij is age in years of the i-th plot observed in the j-th measurement; r is the dimension equal to the number of parameters with fixed effects (global), e_ij is the vector of error terms.

In the fitting of the MEM, the autocorrelation was corrected by modelling the error term (e_ij ) using an autoregressive structure moving average of first order (corARMA1-1). In this study, the dependence between age and DH was modelled in each tree to achieve independence in the residuals of the equation. To verify that corrections for autocorrelation, were appropriate, the likelihood ratio test was performed between model with and without correction, and the adjustment statistics were considered Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC).

The correspondence between the stand information in the polygon mesh and the LiDAR metrics was verified, discarding polygons that did not present recorded returns above 3 m. Finally, 20,098 squares of 400 m² with LiDAR metrics were retained, representing a total of 803.92 ha studied.

Through map algebra, the DH and SI mapping was generated using models L1, M1 and M3 in QGis software version 3.28.11 (QGIS Development Team, 2023). However, some discrepancies were observed; for example, models M1 and M3 presented inconsistencies related to the age assigned to each stand, due to the fact that certain areas were harvested in 2012 under the silvicultural treatment of regeneration method of seed trees. Considering the above, the age of these areas was recalculated by adding to the theoretical age, the age at which harvesting occurred (40 years for P. patula according to the FMP).

To visualize the ranking of the stands based on their potential productivity, they were classified using the most frequent SI. For this, the square polygons with SI within each stand were isolated and the mode was calculated. Subsequently, this value was assigned to each corresponding stand.

The LiDAR-derived metrics with the highest correlation with DH were the 90th, 95th and 99th height percentiles and top height, with similar coefficients (0.96 to 0.97) and linear trend (Figure S1 [suppl.]). The stepwise regression methodology defined p99 as the variable that provided the best fit for predicting DH in the models (Table 3).

Although models L1 and L3 presented outstanding characteristics in terms of fit and similar predictive ability, model L1 (Figure 2) showed superior performance in the normality and homoscedasticity tests of the residuals, according to the Shapiro-Wilk and Breusch-Pagan tests, respectively. Consequently, the L1 model was chosen as the most suitable one to estimate DH based on p99. On the square polygon mesh, the DH prediction was made with the L1 model and p99 and the DH mapping was generated, showing the presence of heights higher than those recorded in the field (35.5 m) (Figure 3).

The DH growth models demonstrated an explanatory power of more than 91% in relation to the total variance of DH growth as a function of age. A margin of error between 1.47 and 1.87 meters was observed. The estimators of the variables were statistically significant at a 95% confidence level (p<0.0001) (Table 4) but the residuals showed autocorrelation, according to the graph ACF.

The goodness-of-fit statistics (RMSE, R adj 2 , AIC and BIC) showed that the M8 model was the most appropriate to model DH. However, the family of growth curves generated (Figure S2 [suppl.] h) overestimates DH with the high SI; this drawback also occurs in model M9 (Figure S2 i) and models M6 (Figure S2 f) and M2 (Figure S2 b). On the contrary, model M5 (Figure S2 e) underestimates the high SI compared to what was observed in the field.

The anamorphic models (M1, M4 and M7) and polymorphic model M3, improved the growth pattern of the SI curve families; however, models M1 (Figure S2 a) and M4 (Figure S2 d) slightly overestimated DH at the high SI. On the other hand, models M3 (Figure S2 c) and M7 (Figure S2 g) corrected this problem; however, both models showed very similar predictions for all SIs at young ages. Consequently, the statistics of these four models suggested their indifferent use; however, when considering the goodness-of-fit statistics and the DH growth curves, the models M1, M3, M4 and M7, were the most appropriate to represent the DH growth pattern of P. patula in the managed forests of SMIC-Atopixco.

The M1-A, M3-A, M4-A and M7-A models, corresponding to the M1, M3, M4 and M7 models respectively, adjusted as global models under the MEM approach with autocorrelation correction, showed high significance in their parameters (Table 5), in addition the SI curves adequately describe the growth pattern in DH (Figure 4). The goodness-of-fit statistics of each model were considered and it was determined that the predictive capacity of the M1-A model is superior compared to the rest of the models. In addition, in the graphical analysis of the standardized versus predicted residuals, an homoscedastic behaviour without abnormal behaviour is observed, while the expected value of the residuals tends to zero (Figure 5 a) and the ACF graphs showed no autocorrelation trend (Figure S3 [suppl.]).

The global models M1-A and M4-A were integrated into each square polygon of the generated mesh to determine the SI. The L1, M1-A and M4-A models can be easily applied to each polygon of the generated mesh. For the estimation of the DH, the L1 model requires that each polygon has the associated p99 and, subsequently, the M1-A and M4-A models require the predicted DH to consider these values as HD₁. For example, with the use of the M1-A model in a polygon that is in a stand of age 25 years, with p99 of 25.63m and base age of 40 years, the projected DH will be given by the expression:

Since the definition of SI is the numerical value given by the DH of a stand at a given age (base age) (Alder, 1980), it is considered that the DH₂ and SI at the age 40 years are equivalent.