Introduction 

Sustainably managed forest ecosystems provide a wide range of social, economic, and environmental benefits, among which carbon dioxide (CO₂) sequestration is particularly significant ([11]). Through photosynthesis, vegetation captures atmospheric CO₂ and converts it into biomass distributed across multiple pools, including living aboveground biomass (stems, branches, foliage), belowground biomass (roots), dead organic matter (litter and coarse woody debris), and soil organic carbon. This continuous carbon accumulation process contributes to climate change mitigation ([15]). In recent years, growing attention has focused on forests’ contribution to biogeochemical cycles, particularly carbon dynamics and their interactions with greenhouse gases ([30]).

Quantifying forest biomass is crucial because it reveals both the carbon sequestration potential of an area and the carbon emissions resulting from its combustion ([28], [8]). In addition, biomass plays a pivotal role in ecosystem productivity and nutrient cycling ([29]) and serves as a valuable ecological indicator of forest structure and sustainability ([4]). It also provides basic information for estimating carbon reservoirs and modeling forest dynamics. Accurate biomass estimation is also a key requirement for Emissions Trading Systems (ETSs), which aim to remove or avoid greenhouse gas emissions ([6]). Under conditions of heterogeneous tree growth, or when new tree species are considered, developing allometric equations from destructive sampling is necessary to obtain accurate biomass stock estimates ([33]).

Carbon storage capacity in biomass depends on stand characteristics such as species composition, age, and density of tree strata ([13]). Although young trees typically exhibit higher rates of carbon assimilation, mature stands store high total carbon stocks, mainly concentrated in the stem ([23], [51]). Conifers allocate a significant proportion of their biomass to the trunk, which represents the main contributor to total tree biomass. For instance, in Abies religiosa (Kunth) Schltdl. & Cham., approximately 84.5% of the biomass is allocated to the stem, which contains about 46% carbon ([1]). Variations of around 45% reported in other studies are attributed to differences in site conditions and methodological approaches ([41], [42], [43]). Similarly, in some Pinus species, aboveground carbon contents ranging from 50% to 53 % have been documented, influenced by stand structural attributes (e.g., basal area and density) as well as edaphic and climatic factors ([44]).

Forest biomass encompasses the weight of organic matter above and below ground and serves as a proxy for ecosystem productivity ([56]). It is commonly partitioned into components such as stem, branches, leaves, bark, and roots ([54], [18]). Tree biomass estimation can be performed through direct or indirect approaches. The direct method involves destructive sampling, which consists of felling trees and weighing their components to obtain precise biomass measurements ([35], [13]). In contrast, indirect methods use allometric equations developed from linear or nonlinear regression analyses using dendrometric variables such as diameter at breast height and total height ([52], [38]), or combinations of variables like diameter at breast height, height, crown cover, leaf area, and wood density ([35]).

Recent studies indicate that species-specific allometric equations reduce bias in above-ground biomass estimates compared with generalized models, which often overlook interspecific and environmental variability and increase uncertainty ([26], [57]). Therefore, it is recommended to use species-specific local equations derived from forest inventory data to improve the accuracy of biomass and carbon estimates ([25]).

Although Vargas-Larreta et al. ([53]) developed national-level allometric equations for 97 forest species from temperate and tropical ecosystems, this dataset still needs to be expanded to include species of high economic and ecological importance. Previous studies on Abies religiosa ([16], [1], [42], [43]) have contributed valuable information, yet their limited scope highlights the necessity to expand the available databases for forestry and silviculture applications.

The genus Abies is the second most economically important in the family Pinaceae in Mexico, contributing 2.8% to the country’s annual timber production ([49]). Its distribution covers only 0.07% to 0.16% of the national territory and is characterized by fragmented and isolated stands ([22]). The present study aims to generate new information on aboveground biomass, carbon content, and carbon dioxide equivalent (CO₂e) of Abies religiosa. This effort responds to both the species’ economic relevance and the current data gaps across Regional Forest Management Units (UMAFORs) in Mexico. Such information is essential to support sustainable management practices, contribute to climate change mitigation, and improve the environmental services provided by Abies forests.

The objective of the study was to develop additive allometric equation systems to accurately estimate aboveground biomass, carbon, and CO₂e across UMAFORs in Mexico. The proposed hypothesis states that biomass production, carbon content, and carbon equivalent do not differ significantly among UMAFORs.

  Materials and methods 

Study area

Six states were considered in this study, covering 13 Regional Forest Management Units (UMAFORs - Fig. 1), distributed as follows: Guerrero (UMAFOR no. 1203), Hidalgo (no. 1303), State of Mexico (no. 1503, 1508, 1509, 1510), Michoacan (no. 1604, 1605, 1608), Puebla (no. 2101, 2105) and Tlaxcala (no. 2901, 2902).

Fig. 1 - Location of the Regional Forest Management Units (UMAFORs) in Central Mexico.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Data collection

A comprehensive sampling of 1919 trees was conducted across managed and unmanaged areas within each UMAFOR, covering the full range of existing diameter classes. In managed areas, we employed a destructive sampling method that involved felling, sectioning, and careful measurement of healthy Abies religiosa trees, while excluding individuals with visible defects that could skew the data. Conversely, non-destructive sampling was carried out in unmanaged areas using advanced laser calipers (Haglof Digitech Professional, Sweden) to measure stem diameters at multiple heights (0.30, 1.30, 2.50 m - [50]).

For each sampled tree, the diameter at breast height (DBH, cm) was measured at 1.30 m using a diameter tape (Forestry Suppliers, Inc., Jackson, MS, USA), while a Haga hypsometer was employed to measure the height (m) of standing trees.

The total length (m) of each felled tree was measured with a Uline Accu-Lock H-1766 tape (WI, USA). Initial measurements were taken from two 0.30 m sections above the stump, followed by an additional measurement at breast height (1.30 m). Subsequent sections were cut at uniform intervals of 2.54 m along the stem, continuing this process until reaching the tree top. For branch biomass and carbon quantification, only branches with a diameter greater than 7.5 cm were included, which corresponds to the minimum merchantable diameter established in the Temperate Forest Planning System (Sistema de Planeación Forestal de Bosque Templado) guidelines ([50]). Tab. S1 (Supplementary material) shows the variables used in the corresponding analyses.

Compilation of information

The stem sections were assumed to have either a cylindrical or paraboloid shape, and their volumes were calculated using the Smalian formula (eqn. 1). For the final section (the tree top), the volume was determined using the cone formula (eqn. 2 - [48]). The overall volume of the stem (V, m³) was considered to be the sum of the volumes of all individual sections (eqn. 1, eqn. 2):

where S₀ is the initial diameter of the section (cm), S₁ is the final diameter of the section (cm), L is the length (m), and V is the volume (m³).

Data analysis and processing

To calculate dry biomass, a basic wood density value of 360 kg m^-3 was used ([16]). From this biomass, the carbon content was estimated, assuming a concentration of 45 % in the dry matter ([41], [42], [43]). To estimate the amount of carbon dioxide equivalent (CO₂e) emitted to the atmosphere from biomass burning, the following formula was applied (eqn. 3):

CO₂e is derived by dividing the molecular weight of CO₂ (44.01 g mol^-1) by the atomic weight of C (12.011 g mol^-1 - [39]).

Systems of additive biomass equations evaluated

The evaluation of additive biomass equations included linear and nonlinear allometric models, following the approaches described by Huy et al. ([20]) and Ordóñez-Prado et al. ([38]). This exploratory analysis aimed to identify the most suitable models and additive equation systems for modeling aboveground biomass and carbon in the UMAFORs included in the study. Based on model performance criteria, two additive equation systems, referred to as S1 and S2 (Tab. 1), were pre-selected. These systems showed the best predictive performance for estimating biomass and aboveground carbon at both the structural component level (stem and branches) and the whole-tree level within each UMAFOR.

Equations fitting method

We applied the dummy variable technique ([32], [14]) to the complete database of each UMAFOR to determine whether specific parameters in the additive equation systems required adjustment. This approach improves model flexibility by incorporating indicator variables that capture additive effects associated with discrete parameters ([38]). We assessed the statistical significance of these parameters at α = 0.05 to determine whether a given UMAFOR required specific parameter values. When no significant differences were detected, a common parameter value was applicable to all UMAFORs ([32]).

To account for potential correlations among the component equations, the Seemingly Unrelated Regression (SUR) method was applied, which allows for the simultaneous estimation of aboveground biomass and carbon components ([40]). Additionally, equation-specific weighting factors were incorporated to correct for heteroscedasticity. We fitted each system of additive equations using the Weighted Nonlinear Seemingly Unrelated Regression (WNSUR) technique, which minimizes residual error, maximizes the statistical significance of the parameter coefficients, and assumes normally distributed and independent errors ([55]).

We fitted each system simultaneously using the MODEL procedure in SAS/ETS v. 9.3 ([47]). After assessing parameter significance using the dummy variable approach, we retained only statistically significant parameters (p < 0.05) in the final models. We corrected heteroscedasticity using weighting factors derived from the error variance structure. Among the tested alternatives, the weighting factor 1/(D²H) (where D is DBH and H is tree height) provided the most effective correction for residual heteroscedasticity ([21], [20]). We performed the Shapiro-Wilk normality test to verify whether the residuals met the normality assumption.

Evaluation and validation of equations

We evaluated the goodness of fit of the additive systems using the adjusted coefficient of determination (R²-adj) and the root mean squared error (RMSE), according to Huy et al. ([20]) and Ordóñez-Prado et al. ([38]). To compare the predictive performance of the S1 and S2 systems, we applied the “k-fold” cross-validation method ([24]), which allows for robust model comparison and selection of the best-performing system ([20]). For cross-validation, we calculated error metrics such as Bias (%), root mean square error (RMSE, %), and mean absolute percentage error (MAPE, %). We calculated these statistics ten times and averaged the results. We also examined diagnostic plots of predicted vs. observed values and weighted residuals vs. fitted values to evaluate the equation’s behavior. We considered the equations with the lowest cross-validation error values the most suitable ([21] - eqn. 4, eqn. 5, eqn. 6):

where k is the number of iterations (k=10), n is the number of trees sampled for validation, y_i and hat{y}_i are biomass or carbon per component and total observed and predicted for the i-th tree sampled at iteration k, respectively. The final parameter estimates for the selected equation systems were obtained by fitting the entire dataset using the MODEL procedure in SAS/ETS v. 9.3 ([47]).

  Results 

The estimation of aboveground biomass and carbon revealed a clear structural differentiation in the contribution of individual tree components. The stem (Bst) was identified as the primary reservoir (97%), representing the largest share of total aboveground biomass (AGB). In contrast, branches (Bbr) accounted for a smaller fraction (3%) that was still significant in the overall accumulated carbon stock (AGC).

Significant contrasts were observed among the evaluated UMAFORs (Tab. 3), in terms of extreme differences (minimum and maximum values) of biomass components (Bst, Bbr, and AGB) and their carbon analogs (Cst, Cbr, and AGC). This marked inter-UMAFOR variability directly reflects the inherent heterogeneity in the composition and structure of the stands evaluated (see Tab. S1 in Supplementary material).

Goodness-of-fit and selection of the best system

The results showed that system S2 achieved higher accuracy in estimating the evaluated variables compared with system S1, which exhibited higher RMSE and RMSE% values. Although both systems yielded high R²-adj values, S2 was slightly superior, indicating a more efficient explanation of data variability. Moreover, S2 had lower Bias% and MAPE% values, indicating a reduced tendency to under- or overestimate predictions. Consequently, S2 yielded more consistent predictions and a closer alignment with observed values, offering higher predictive reliability for both biomass components and total aboveground carbon (Tab. 2).

The specific and statistically significant parameters (p<0.05) obtained for each UMAFOR through the fitting of the additive equations using the WNSUR procedure and the dummy variable technique for the system S2 are presented in Tab. S2 (Supplementary material). System S2 stands out as the most reliable and robust option, as evidenced by its low standard error values. All estimated parameters were statistically significant (p<0.05), which confirms the model suitability for estimating aboveground biomass and carbon at both the component level (stem and branches) and the total trees level, as well as for assessing total biomass and carbon across UMAFORs in each state.

Developed systems of equations

The residuals for biomass estimation in both systems were below 0.50 kg, indicating that the equations (S1 and S2) effectively captured biomass variation in the analyzed components. System 2 showed a superior ability to capture the relationship between predicted and observed values, leading to more accurate predictions. The enhanced performance of S2 may be due to better model parameterization, the inclusion of more informative predictor variables, and a higher capacity to represent the intrinsic complexity of biomass distribution. The enhanced predictive performance of S2 is evident from the closer alignment of its regression line to the 1:1 identity line in the total biomass plots (Fig. 2II.f, Fig. 3II.b).

Fig. 2 - Residual analysis and comparison of predicted vs. observed values for aboveground dry biomass (kg tree^-1) by component and total.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Fig. 3 - Residual analysis and comparison of predicted vs. observed values for aboveground carbon (kg tree^-1) by component and total.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Dispersion increased slightly for the branch component (Fig. 2.II.c) in the S2 system, although it remained the best-performing system overall. Estimating branch biomass is inherently challenging, so higher errors were anticipated, even with the best equation. Although the S2 system may exhibit greater errors in branch biomass estimation, it outperformed in accurately estimating total tree aboveground biomass (Fig. 2). Similarly, the graphs depicting carbon estimation confirmed S2 as the most accurate model (Fig. 3II.c), as reflected by the stronger overlap between the observed and predicted value lines across the datasets of both models (Fig. 3).

Estimations at UMAFOR level

The assessment of total aboveground biomass in each UMAFOR revealed significant variability in carbon storage capacity, primarily influenced by differences in tree diameter and stand density within each forest region. The UMAFORs with the highest mean biomass per tree were Puebla-2105 (4252 kg), State of Mexico-1503 (3970 kg), Michoacan-1608 (3783 kg), and Puebla-2101 (3928 kg), with some individual trees storing over 1900 kg of carbon in the most extreme cases. These UMAFORs are characterized by maximum DBH exceeding 80 cm and a mean DBH of 40 cm or more, which contributes to substantial biomass accumulation per individual (Tab. S1 in Supplementary material).

Conversely, the UMAFORs in Guerrero (1203 with 2460 kg) and Tlaxcala (2901 with 1823 kg and 2902 with 1908 kg) exhibited the lowest values for both biomass and DBH, reflecting a lower carbon sequestration capacity at the individual tree level. Nevertheless, some of these UMAFORs, such as 2901 (820 kg) with 202 individuals, still have a considerable number of trees, which may enhance their total carbon contribution despite the lower biomass per tree (Tab. S1). This suggests that, although tree number is a determinant factor of total carbon storage at the management unit level, it plays a comparatively minor role in explaining carbon accumulation per individual tree.

Overall, UMAFORs with larger DBH, such as State of Mexico-1508 (1519 kg), Puebla-2101 (1768 kg), and Michoacan-1608 (1703 kg), and a high number of individuals demonstrated a greater potential for carbon sequestration at the landscape level. While the direct positive relationship between DBH and biomass was evident, tree density also played an essential role in estimating total carbon storage within a single UMAFOR (Tab. 3; see also Tab. S1 in Supplementary material).

Total aboveground biomass varied considerably among forest regions, reflecting differences in carbon storage capacity and potential CO₂e emissions to the atmosphere. This variability is particularly relevant when assessing each UMAFOR’s contribution to carbon storage and release. In terms of CO₂e content, the UMAFOR of Puebla (2105), which has a maximum biomass value of 4252 kg per tree, could emit up to 7803 Mg CO₂e per individual in cases of biomass loss or combustion. Other UMAFORs, such as Puebla-2101, State of Mexico-1503, State of Mexico-1508, and Michoacan-1608, have trees with high aboveground biomass ranging from 3374 to 3970 Mg, corresponding to potential emissions of 6.20 to 7.30 Mg CO₂e per tree under loss scenarios. In contrast, UMAFORs of Tlaxcala (2901 and 2902), with lower maximum biomass values (1822 and 1908 kg), would release approximately 3.3 to 3.5 Mg CO₂e per tree (Tab. S1).

When these values are scaled to the number of trees per hectare under different stand densities (low: 100 trees ha^-1; medium: 250 trees ha^-1; high: 400 trees ha^-1), the potential CO₂e emissions per unit area become even more evident. For example, the UMAFOR Puebla-2105, under medium density, could emit up to 1950.85 Mg CO₂e ha^-1 (Tab. 3). Other UMAFORs, such as State of Mexico-1503 and Puebla-2101, also exceed 1800 Mg ha^-1. Even in low-density areas, these regions maintained a high emission potential, exceeding 700 Mg CO₂e ha^-1. In comparison, the UMAFORs Tlaxcala-2901 and Tlaxcala-2902, with lower aboveground biomass per tree, projected emissions below 875 Mg CO₂e ha^-1 under medium-density conditions (Tab. 3).

  Discussion 

The fitting strategy based on the Weighted Nonlinear Seemingly Unrelated Regression (WNSUR) method, which combines heteroscedasticity correction with the dummy variable technique, demonstrated remarkable statistical robustness in generating additive equation systems at the level of Regional Forest Management Units (UMAFOR) in Mexico. The reliability and predictive capability of the models were validated through a rigorous k-fold cross-validation procedure, which confirmed their stability and low prediction error variance. Integrating heteroscedasticity correction within the WNSUR provided robust standard errors and consistent, asymptotically normal estimators ([40], [21]).

This methodology has been successfully applied in several contexts ([37], [38], [2]), showing high accuracy, reliability, and flexibility in constructing equations consistent with the additivity principle. This principle, fundamental in biomass and carbon modeling, ensures that the sum of estimates for individual tree components (stem, branches, bark, and foliage) is consistent with total biomass ([9]). Recent studies in tropical and temperate forests indicate that this approach improves predictive performance compared with less complex or non-additive models ([20]).

The S2 system showed slightly higher precision and consistency in estimating biomass and aboveground carbon than the S1 system, as shown by lower RMSE, RMSE%, Bias%, and MAPE% values. These metrics suggest reduced error dispersion and minimal systematic bias in S2. Moreover, the applied fitting strategy effectively stabilized residual variance, preventing increases in errors associated with tree size ([40], [10]). Although both systems achieved a high adjusted coefficient of determination (R²-adj > 0.99), S2 more effectively accounted for variability in total aboveground biomass and carbon data at a highly precise scale, corroborating the findings by Hernández-Moreno et al. ([19]) in Michoacán, Mexico.

The slight variability observed in the branch component estimates does not undermine the overall efficacy of the S2 system (Fig. 2.II.c, Fig. 3.II.c), as branch biomass is inherently more variable due to environmental and biological factors ([17]). This behavior is expected, given the structural irregularity of branches, which increases the complexity of biomass estimation. Nevertheless, the S2 system achieved substantially higher accuracy in predicting total aboveground biomass. This assertion is further supported by the k-fold cross-validation results, which confirmed the model’s predictive generalization to independent datasets, a critical metric of model robustness ([27], [45]).

Incorporating total tree height alongside DBH significantly enhanced the accuracy of biomass estimates. This finding highlights the role of height in tree allometry, as trees with similar diameters but different taper forms exhibit different allometric relationships ([52]). Moreover, height is a key indicator of inter-tree competition and site quality, providing additional information that improves the predictive performance of biomass equations ([46]).

Additivity in biomass and carbon equations is essential for ensuring internal consistency between individual tree components and their total estimates ([9]). Vargas-Larreta et al. ([54]) and Flores-Medina et al. ([12]) highlighted the importance of this property, as it prevents inconsistencies between the sum of component estimates and the corresponding total values. However, the application of additive modeling approaches in Mexico remains limited. Notable contributions include Monroy-Rivera & Návar-Cháidez ([31]), who developed additive biomass equations for Hevea brasiliensis (Willd. ex A.Juss.) Müll.Arg. plantations in Veracruz, and Návar et al. ([34]), who estimated biomass for Pinus durangensis Martínez, P. cooperi C.E. Blanco, and P. engelmannii Carr. in Durango, Mexico. In addition, Flores-Medina et al. ([12]) reported satisfactory results for four naturally regenerating species, while Vargas-Larreta et al. ([54]) achieved high model performance for 17 forest species in temperate regions of northwestern Mexico. More recently, Cuevas-Cruz & Aquino-Ramírez ([8]) demonstrated that DBH and total height explained over 92% of the biomass variability across components and 98% for total biomass.

Our findings are consistent with previous results reported in the literature. For example, Avendaño-Hernández et al. ([1]) reported that Abies religiosa in Tlaxcala allocates approximately 84.5% of its biomass to the stem and 6.9% to the branches, with a carbon content of 46.48%. These results are consistent with the high coefficients of determination (R² = 99%) reported by Razo-Zárate et al. ([42]) for biomass and carbon estimations in post-fire stands of Hidalgo. A balanced stem-to-branch biomass ratio was observed within the UMAFORs Tlaxcala-2901, Tlaxcala-2902, and Hidalgo-1303, where comparable R² values (0.99 and 0.97, respectively) were achieved despite the use of different estimating equations.

In UMAFOR 1303, significant contributions to biomass, carbon, and CO₂e were recorded, totaling approximately 493.6 Mg of biomass, 246.9 Mg of carbon, and 906.8 Mg of CO₂e for the sampled 226 trees (Tab. S1 in Supplementary material). These results are comparable to the carbon sequestration values reported by Razo-Zárate et al. ([42]) in fire-affected areas, where 166.68 Mg C ha^-1 (611.2 Mg CO₂e ha^-1) were estimated in disturbed stands and 62.68 Mg C ha^-1 (229.7 Mg CO₂e ha^-1) in conserved sites.

These discrepancies may reflect differences in growth conditions, tree age, and sampling methods. According to Návar-Cháidez & Domínguez-Calleros ([36]) and Cuevas-Cruz et al. ([7]), age- and climate-driven growth rates affect basic wood density, which helps explain these variations. Studies on Abies religiosa in Hidalgo, Mexico, by Chávez et al. ([5]) further indicate that carbon fixation through photosynthesis is higher in younger stands than in mature ones. Moreover, UMAFORs of the State of Mexico demonstrated outstanding carbon sequestration capacities, with CO₂e values ranging from 1318 to 1938 Mg (Tab. S1 in Supplementary material). In contrast, other areas, such as Texcoco, Mexico, showed significantly lower carbon concentration, with an average carbon density of 376 Mg ha^-1 across five evaluated sites, where aboveground biomass is the primary reservoir and accounts for 59% of total carbon stocks ([3]). Although the forested areas of Texcoco show moderate carbon sequestration potential, mature Abies religiosa forests have a notable advantage because of their higher carbon concentration in aboveground biomass.

  Conclusions 

The allometric equations developed in this study provide a robust framework for estimating aboveground biomass, carbon, and carbon dioxide equivalent at the UMAFOR level in Mexico. The slightly superior performance of system S2 highlights its statistical robustness and suitability for national-scale applications involving Abies religiosa. The results emphasize the significant influence of stand structural attributes, such as tree DBH and density, on carbon storage potential. The additive property of the proposed equations ensures consistency between individual component estimates and total tree biomass, facilitating their practical application by forestry professionals. Consequently, these allometric equations represent valuable tools for forest monitoring, sustainable management, and accurate quantification of carbon stocks. Furthermore, they support climate change mitigation strategies aimed at conserving and managing Abies religiosa forests across Mexico.

  Acknowledgements 

The first author gratefully acknowledges the Secretaría de Ciencias, Humanidades, Tecnología e Innovación (SECIHTI) for providing the postdoctoral fellowship at the Colegio de Postgraduados, Campus Montecillo, Texcoco. We also thank the agrarian authorities of the community of San Juan Quiahije, Juquila, and the forestry technicians, Gregorio Morales and María del Sol, for their assistance and support during the fieldwork.

  References

  Supplementary Material

Tab. S1 - Descriptive statistics of the data used.
Tab. S2 - Specific parameters by UMAFOR obtained from the fitting of the additive equation system with the WNSUR procedure and the dummy variable technique using system S2.