Study area

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Ángeles-Pérez G, Méndez-López B, Valdez-Lazalde JR, Plascencia-Escalante FO, de los Santos-Posadas HM., Chávez-Aguilar G, Ortiz-Reyes , AD, Soriano-Luna , MÁ, Zaragoza-Castañeda Z, Ventura-Palomeque E, Martínez-López A, WaysonC, López-Merlín D, Olguín-Álvarez M, Carrillo-Negrete O, Maldonado-Montero V, 2015. Estudio de caso del Sitio de Monitoreo Intensivo del Carbono en Hidalgo. Colegio de Postgraduados, Texcoco, México.) 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., 2021Pérez-Vázquez ZR, Ángeles-Pérez G, Chávez-Vergara B, Valdez-Lazalde JR, Ramírez-Guzmán ME, 2021. Enfoque espacial para modelación de carbono en el mantillo de bosques bajo manejo forestal maderable. Madera Bosques 27(1). 10.21829/myb.2021.2712122), 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., 2015Ortiz-Reyes AD, Valdez-Lazalde JR, de los Santos-Posadas HM, Ángeles-Pérez G, Paz-Pellat , F, Martínez-Trinidad T, 2015. Inventario y cartografía de variables del bosque con datos derivados de LiDAR: comparación de métodos. Madera Bosques 21(3): 111-128. 10.21829/myb.2015.213461).

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., 2010Cruz-Leyva IA, Valdez-Lazalde JR, Ángeles-Pérez G, de los Santos-Posadas HM, 2010. Modelación espacial de área basal y volumen de madera en bosques manejados de Pinus patula y P. teocote en el ejido Atopixco, Hidalgo. Madera Bosques 16(3): 75-97. 10.21829/myb.2010.1631168). 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, 2004Cruz-Ruiz F, 2004. Programa de Manejo Forestal para el aprovechamiento de recursos forestales maderables del Ejido Atopixco. Zacualtipán de Ángeles, Hidalgo, México.).

Light Detection and Ranging-derived data

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., 2015Ortiz-Reyes AD, Valdez-Lazalde JR, de los Santos-Posadas HM, Ángeles-Pérez G, Paz-Pellat , F, Martínez-Trinidad T, 2015. Inventario y cartografía de variables del bosque con datos derivados de LiDAR: comparación de métodos. Madera Bosques 21(3): 111-128. 10.21829/myb.2015.213461; Galeote-Leyva et al., 2022Galeote-Leyva B, Valdez-Lazald e JR, Ángeles-Pérez G, de los Santos-Posadas HM, Romero-Padilla JM, 2022. Inventario asistido por LiDAR: efecto de la densidad de retornos y el diseño de muestreo sobre la precisión. Madera Bosques 28(2). 10.21829/myb.2022.2822330). 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, 2023R development Core Team. (2023). R: A Language and Environment for Statistical Computing. (4.3.1). R Foundation for Statistical Computing, Vienna, Austria.). 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, 2023QGIS Development Team, 2023. QGIS geographic information system (3.28.11). Open-Source Geospatial Foundation.) 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.

Dendrometric data

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 (2008Hollinger DY, 2008. Defining a Landscape-Scale Monitoring Tier for the North American Carbon Program. In: Field Measurements for Forest Carbon Monitoring. A Landscape-Scale Approach; HooverCM (ed). pp: 3-16. Springer, Netherlands. 10.1007/978-1-4020-8506-2_1). Site data were processed assuming spatial independence among the study variables.

  

Figure 1 Components of the Intensive Carbon Monitoring Site (SMIC) Atopixco, a) Mexico, b) Zacualtipán de Ángeles, and c) study area and sampling scheme. 

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 (1999Bengoa-Martínez JL, 1999. Estimación de la altura dominante de la masa a partir de la “altura dominante de parcela”: ventajas frente a la altura dominante de Assman. For Syst 8(3): 311-321.)), which were incorporated into the database with the LiDAR metrics per site.

Fitted models for site index

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., 1983Clutter JL, FortsonJC, Piennar LV, Brister GH, Bailey RL, 1983. Timber management: a quantitative approach. John Wiley & Sons, Inc.; Cieszewski & Bailey, 2000Cieszewski CJ,Bailey RL, 2000. Generalized Algebraic Difference Approach: Theory Based Derivation of Dynamic Site Equations with Polymorphism and Variable Asymptotes Dynamic Equation Modeling Background. Forest Sci 46(1): 116-126. 10.1093/forestscience/46.1.116) (Table 2). Each of the models was structured in the form DH₂, $=f\left({\text{E}}_{\text{1}}\text{,}{\text{E}}_{\text{2}}\text{,}{\text{DH}}_{\text{1}}{\text{β}}_{\text{i}}\right)$, where ${\text{DH}}_1$ is the DH to be predicted at age ${\text{E}}_{\text{2}}$; ${\text{β}}_{\text{i}}$ are the parameters to be estimated ${\text{DH}}_{\text{1}}$ and ${\text{E}}_{\text{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., 1983Clutter JL, FortsonJC, Piennar LV, Brister GH, Bailey RL, 1983. Timber management: a quantitative approach. John Wiley & Sons, Inc.). 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, 2023R development Core Team. (2023). R: A Language and Environment for Statistical Computing. (4.3.1). R Foundation for Statistical Computing, Vienna, Austria.) 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 $\left({\text{R}}_{\text{adj}}^{\text{2}}\right)$, 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, 2003Kozak A, Kozak R, 2003. Does cross validation provide additional information in the evaluation of regression models?Can J Forest Res 33(6): 976-987. 10.1139/x03-022), 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 (2001Fang Z, Bailey RL, 2001. Nonlinear Mixed Effects Modeling for Slash Pine Dominant Height Growth Following Intensive Silvicultural Treatments. Forest Sci 47(3): 287-300. 10.1093/forestscience/47.3.287), for this, the following model was created: ${HD}_2=f\left(A_i\beta +{\beta }_i{b}_i,t_{ij}\right)+{\epsilon }_{ij}$, 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).

Site index mapping

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, 2023QGIS Development Team, 2023. QGIS geographic information system (3.28.11). Open-Source Geospatial Foundation.). 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.

Dominant height models with LiDAR

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).

  

Figure 2 Graphical fit of the L1 model a) regression line, b) standardized residuals in each sampling subunits in the Intensive Carbon Monitoring Site (SMIC) Atopixco. 

  

Figure 3 Mapping of the dominant height (DH) estimated by the L1 model per pixel at the Intensive Intensive Carbon Monitoring Site (SMIC) Atopixco. 

Families of SI curves

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, ${\text{R}}_{\text{adj}}^{\text{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.]).

  

Figure 4 Site index curves for Pinus patula at the base age of 40 years under the mixed-effects method (MEM), a) M1-A model, b) M3-A model, c) M4-A model and d) M7-A model in the Intensive Carbon Monitoring Site (SMIC) Atopixco. 

  

Figure 5 Standardized residuals vs. predicted from global models a) M1-A model, b) M3-A model, c) M4-A model and d) M7-A model in the Intensive Carbon Monitoring Site (SMIC) Atopixco. 

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, 1980Alder D (ed), 1980. Forest Volume Estimation and Yield Prediction: Yield prediction: Vol. FAO Forestry Paper (2nd ed.). Food and Agriculture Organization of the United Nations, Rome. 194 pp.), it is considered that the DH₂ and SI at the age 40 years are equivalent.