Allometric Equations for Estimating Biomass and Carbon Stocks in the Temperate Forests of North-Western Mexico

This paper presents new equations for estimating above-ground biomass (AGB) and biomass components of seventeen forest species in the temperate forests of northwestern Mexico. A data set corresponding to 1336 destructively sampled oak and pine trees was used to fit the models. The generalized method of moments was used to simultaneously fit systems of equations for biomass components and AGB, to ensure additivity. In addition, the carbon content of each tree component was calculated by the dry combustion method, in a TOC analyser. The results of cross-validation indicated that the fitted equations accounted for on average 91%, 82%, 83% and 76% of the observed variance in stem wood and stem bark, branch and foliage biomass, respectively, whereas the total AGB equations explained on average 93% of the total observed variance in AGB. The inclusion of total height (h) or diameter at breast height2 × total height (d2h) as a predictor in the d-only based equations systems slightly improved estimates for stem wood, stem bark and total above-ground biomass, and greatly improved the estimates produced by the branch and foliage biomass equations. The predictive power of the proposed equations is higher than that of existing models for the study area. The fitted equations were used to estimate stand level AGB stocks from data on growing stock in 429 permanent sampling plots. Three machine-learning techniques were used to model the estimated stand level AGB and carbon contents; the selected models were used to map the AGB and carbon distributions in the study area, for which mean values of respectively 129.84 Mg ha−1 and 63.80 Mg ha−1 were obtained.


Introduction
Improved knowledge of carbon stocks and fluxes is needed in order to understand the current state of the carbon cycle and how it might evolve with changing land use and climatic conditions [1].This has led to an increased interest in estimating forest biomass for practical forestry purposes and for scientific purposes.Tree biomass is an important component of the carbon pool in forests, and it can be estimated by using biomass expansion factors [2] or by relating biomass functions to tree-level data Forests 2017, 8, 269 2 of 20 obtained in forest inventories [3].In both cases, generic biomass functions are used to quantify the carbon in forests [4] because such functions improve the accuracy of carbon accounting systems and thus allow accurate planning of whole-tree and residual biomass utilization for bioenergy production.
Mexico is one of the largest emitters of CO 2 due to deforestation and land degradation, which together account for 4.9 per cent (32 Mt CO 2 e) of the global emissions in the country [5].Thus, most states in Mexico are implementing action plans aimed at mitigating the effects of climate change and accessing economic incentives that favour carbon sequestration in forests.Initiatives such as REDD+ (Reducing Emissions from Deforestation and Forest Degradation and enhancement of carbon stocks) are important efforts aimed at combating climate change; however, for effective implementation of such mechanisms, accurate estimation and monitoring of AGB and associated carbon stocks in forests is first required [6].
Several biomass equations and biomass expansion factors have been reported for forest species in Mexico.Acosta-Mireles et al. [7] generated allometric equations for six species of a mixed hardwood and oak forest in the northern Sierra of Oaxaca.All equations were of the form Y = bX k , where Y is the above-ground biomass (kg), X is diameter at breast height (cm), and b and k are the parameter estimates.Návar [8,9] developed diameter-based allometric biomass equations for nine species of pine and oak in temperate forests in Durango and Chihuahua.Rodríguez-Ortíz et al. [10] analyzed the effects of thinning on distribution and accumulated above-ground biomass content of Pinus patula in Ixtlán, Oaxaca.Silva-Arredondo and Návar-Cháidez [11] conducted research in northern Mexico to assess biomass components and biomass expansion factors by using classical allometric equations, and Villegas-Jiménez et al. [12] estimated above-ground biomass components in 13 Mexican provenances of P. greggii Engelm. in southern Mexico.Other similar studies have been carried out in Mexico by Díaz-Franco et al. [13], Mendoza-Ponce and Galicia [14], Búrquez et al. [15], Aguirre-Salado et al. [16], and Méndez-González et al. [17].Rojas-García et al. [18] provide a comprehensive database with the above cited and other allometric equations applicable to a large number of tree and shrub species in Mexico.However, biomass equations were only available for nine tree species in the study area [8,9].Moreover, these equations were developed using tree diameter as the only regressor, and they are not additive, one of the most important properties of biomass equations.Thus, there is still a need for research involving the accurate estimation of forest biomass and carbon stocks in these multi-species forests.
Temperate forests occupy 32,330,508 ha in Mexico, i.e., 17 per cent of the total land area.These are the richest ecosystems in Mexico with some 7000 plant species [19], of which ~150 are species of pines and 170 of oaks together representing more than 50 per cent of all known pine and oak species [20].Durango is the most important forestry state in Mexico, with 10.5 million ha of forest cover and a timber production of 1.9 million m 3 per year (32.8% of the national forest production [21]).Quantification of forest biomass and carbon stocks is important for management of these forests.Furthermore, mapping the spatial pattern of large-scale forest biomass can provide an overall picture of the carbon stocks within a region, which is of great scientific and political importance [22].However, as pointed out in studies by Zianis and Mencuccini [23] and Walker et al. [24], although forest biomass can be measured directly, the method used (destructive sampling) is expensive and time-consuming and does not yield the continuous spatial distribution of biomass at large scales [25].The difficulties associated with the direct estimation of biomass must be solved by quantitative satellite and remote sensing techniques, which represent the most promising methods for estimating biophysical features of forests, along with the development of parametric and nonparametric statistical methods for modelling stand variables [26].Together these methods have enabled estimation of biomass at multiple scales with large spatial and temporal coverage [27,28].
The main objective of the present study was to develop species-specific systems of additive equations for predicting total above-ground biomass in the temperate forests of north-western Mexico.In addition, we used the new equations to model the forest biomass and carbon distribution at regional Forests 2017, 8, 269 3 of 20 scale by applying remote sensing Landsat-5 TM imagery, terrain parameters and data from permanent research plots.

Study Area
The study was conducted in the uneven-aged and multi-species temperate forest of Durango ( 22• 20 49 to 26 • 46 33 N; 103 • 46 38 to 107 • 11 36 W), which occupies about 23% of the Sierra Madre Occidental ecosystem (Figure 1).The study area is divided in Regional Forest Management Units (in Spanish, Unidades de Manejo Forestal Regional: UMAFORES).The elevation above sea level varies between 363 and 3200 m (average 2264 m).Precipitation ranges from 443 to 1452 mm, with an annual average of 917 mm, whereas the mean annual temperature varies from 8.2 to 26.2 • C, with an annual average of 13.3 • C. The predominant forest types are pine and uneven-aged pine-oak, with one or two pine species (usually P. cooperi and P. durangensis) dominating the overstorey and Q. sideroxyla dominating the understorey [29].
Forests 2017, 8, 269 3 of 20 at regional scale by applying remote sensing Landsat-5 TM imagery, terrain parameters and data from permanent research plots.

Study Area
The study was conducted in the uneven-aged and multi-species temperate forest of Durango (22°20′49″ to 26°46′33″ N; 103°46′38″ to 107°11′36″ W), which occupies about 23% of the Sierra Madre Occidental ecosystem (Figure 1).The study area is divided in Regional Forest Management Units (in Spanish, Unidades de Manejo Forestal Regional: UMAFORES).The elevation above sea level varies between 363 and 3200 m (average 2264 m).Precipitation ranges from 443 to 1452 mm, with an annual average of 917 mm, whereas the mean annual temperature varies from 8.2 to 26.2 °C, with an annual average of 13.3 °C.The predominant forest types are pine and uneven-aged pine-oak, with one or two pine species (usually P. cooperi and P. durangensis) dominating the overstorey and Q. sideroxyla dominating the understorey [29].

Biomass Data
Biomass estimates for stem wood, stem bark, branches and foliage were obtained from 1336 destructively sampled trees of 17 species: Pinus cooperi, P. durangensis, P. engelmannii, P. leiophylla, P. herrerae, P. teocote, P. lumholtzii, P. strobiformis, P. oocarpa, P. douglasiana, P. michoacana, Juniperus depeanna, Arbutus bicolor, Quercus sideroxyla, Q. rugosa, Q. durifolia and Q. crassifolia.As the number of specimens of P. michoacana and P. douglasiana was small, and as these species have similar morphological characteristics and are located in the same forest stands [30], the data for these species

Biomass Data
Biomass estimates for stem wood, stem bark, branches and foliage were obtained from 1336 destructively sampled trees of 17 species: Pinus cooperi, P. durangensis, P. engelmannii, P. leiophylla, P. herrerae, P. teocote, P. lumholtzii, P. strobiformis, P. oocarpa, P. douglasiana, P. michoacana, Juniperus depeanna, Arbutus bicolor, Quercus sideroxyla, Q. rugosa, Q. durifolia and Q. crassifolia.As the number of specimens of P. michoacana and P. douglasiana was small, and as these species have similar morphological characteristics and are located in the same forest stands [30], the data for these species were combined Forests 2017, 8, 269 4 of 20 in order to yield sufficient data for model fitting (Table 1).A generic equation was also developed for each biomass component for all pines and all oak trees.
For above-ground biomass measurement, the fieldwork included selection, measurement and felling of standing trees, collection of dimensional data, cutting and separating the tree components and weighing fresh components of each biomass component on site.The sample trees were subjectively selected to represent the complete range of diameter at breast height and total tree height; thus, the data were representative of the different growth and structure conditions in the study area.The number of trees sampled varied from a minimum of 30 in the case of P. douglasiana to a maximum of 130 trees of P. durangensis and Q. crassifolia.Diameter at breast height (d) and total height (h) were measured in each sample tree.The number of trees sampled per species depended on the distribution and abundance within the study area.
The following biomass components were considered: stem wood, stem bark, branches (including both wood and bark) and foliage (needles/leaves).Stems were divided into smaller sections to facilitate the weighing process.For each felled tree, stem diameter outside bark was measured at 0.3, 0.6, 1.3 m and thereafter every 2.5 m along the stem to take into account variations in moisture content.The total height of each felled tree was measured with a long fiberglass tape.The green weight of stem and branches was determined by weighing the logs and branches in the field by placing them on a 1000 kg balance (precision 100 g).Whole foliage was totally separated from the trunk and weighed on an analytical balance (precision 1 g).A disc of about 5 cm thick was cut from each log, and representative samples of branches and needles/leaves were weighed in the field (fresh weight) before being transported to the laboratory, where they were oven-dried at 70-85 • C to constant weight (dry weight, measured to the nearest 0.1 g).The time required to reach constant weight was 7-9, 3-4 and 2-3 days for wood discs, branches and foliage, respectively.The biomass of each tree component was calculated on the basis of the ratio of dry biomass to fresh biomass and then summed to produce the total AGB of each tree sampled.The carbon content of each tree component was determined by the dry combustion method, in a TOC analyser (Ol Analytical, College Station, Texas, USA).The number of observations used for tree biomass estimations and the basic description of the tree biomass components data for each species are summarized in Table 1.where n = number of trees, d = diameter at breast height (cm), h = total height (m), W w = stem wood biomass (kg tree −1 ), W b = stem bark biomass (kg tree −1 ), W br = wood plus bark biomass of branches (kg tree −1 ), W f = foliage (leaves/needles) biomass (kg tree −1 ), W t = total above-ground biomass (W w + W b + W br + W f ) (kg tree −1 ).

Basic Models
We used three basic model forms as starting points for model selection (Equations ( 1)-( 3)).
where α, β, and γ are the equation parameters, w can be total tree AGB or any of the tree biomass components considered in the study, and ε i is the model error.
In a first step the models were linearized by using natural logarithms, and the stepwise variable selection approach was used to select the best set of independent variables for each biomass component equation.The significance level for entering and maintaining variables in the model was restricted to 0.001 [2].In a second step, a species-specific system of equations with cross-equation constraints on the structural parameters and cross-equation error correlation was then defined for additive prediction of tree component and above-ground biomass [31,32].
The species-specific biomass equation system is formulated as follows: where W i represents the tree biomass for the i-th component, W t is the total tree AGB (i.e., the sum of all the tree biomass components), X j are tree variables (j = d, h, d 2 h), α i , β ij are parameters to be estimated in the fitting process, and ε i , ε t are inter-correlated error terms [32].

Simultaneous Fitting of Tree Biomass Components and Total AGB
Total tree AGB was formulated as the sum of the equations for each tree component, and the system of equations was fitted using the Generalized Method of Moments (GMM) in the PROC MODEL procedure of SAS/ETS® [33].This method produces efficient parameter estimates under potentially heteroscedastic conditions, without specifying the nature of the heteroscedasticity [34], and it thus avoids estimating the heteroscedastic error variance.Bi et al. [32] and Castedo et al. [2] state that the GMM procedure of SAS overcomes the above-mentioned problem of the linear combination of the error terms of the tree biomass component equations by computing a generalized inverse of the variance covariance matrix (by setting part of the matrix for whole-tree biomass to zero) and that estimation of parameters by simultaneously fitting the biomass component equations guarantees that AGB will be the sum of the tree component estimates.

Comparison of Equations
Statistical and graphical analyses were used to compare the performance of the equations.Moreover, a 10-fold cross-validation approach was used to evaluate the predictive performance of the models.The root mean square error (RMSE) and the coefficient of determination (R 2 ), calculated by Equations ( 6) and (7), were estimated using the residuals from the 10 fold cross-validation.where Y ij and Ŷij are the j-th observed and predicted values of biomass for component i, Y ij is the mean of n observed values for the same component, and p is the number of parameters in the model.

Applying the Proposed Above-Ground Biomass Equations
A data set from a network of 429 permanent research sampling plots (Sitios Permanentes de Investigación Forestal y de Suelos: SPIFyS), distributed across the Sierra Madre Occidental in Durango was used to relate stand biomass and carbon stocks to variables obtained from remote sensors.The plots were installed between 2007 and 2011 using the protocol developed by [35].The data from these permanent sample plots are used to monitor the growth and yield of the main types of forest in Durango.The plots are 50 × 50 m in size and they were distributed by systematic sampling, with a variable grid ranging from 3 to 5 km, depending on the size of the forest properties.The sampling plots are intended to be re-measured at five-year intervals.Among other variables, tag number, species code, breast height diameter (measured in cm at 1.3 m above ground level), total tree height (m), height to the live crown (m), azimuth (•) and radius (m) from the centre of the plot to all trees equal or larger than 7.5 cm in diameter were recorded.The total AGB in each stand was calculated by applying the new tree-level biomass equations (to each tree), converted into carbon content by using the carbon proportion estimated for each component and expressed per unit area (ha).The species or group species-specific AGB models reported by Rojas-García et al. [18] were used to estimate total AGB for the tree species present in the permanent plots and for which no biomass equations were developed in this study.
The spectral data were derived from Landsat TM5 (Thematic Mapper) satellite images obtained during the same period that the SPIFyS were established (2007 to 2011) and are freely available from the National Landsat Archive Processing System (NLAPS).The Landsat 4-5 Thematic Mapper product, level 1 of surface reflectance (radiometrically and atmospherically corrected) was processed using the Standard Landsat Product Generation System (LPGS) via the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm (both available at [36]), which uses a 6S radiative transfer model to correct for atmospheric effects on a given date [37].Landsat TM5 bands 1, 2, 3, 4, 5 and 7 were used.Band 6, which is designed for thermal mapping and determining soil moisture, was not considered because of its lower (120 m) spatial resolution [38].Some common Vegetation Indexes (VI) and other derived parameters were computed from the atmospherically corrected image bands: the Normalized Difference Vegetation Index (NDVI); the Soil Adjusted Vegetation Index (SAVI) and its modification (MSAVI); the Enhanced Vegetation Index (EVI); the Normalized Burn Ratio using bands 4 and 7 (NBR1); the difference Normalized Burn Ratio 2 (NBR2) using bands 5 and 7; and the Normalized Difference Moisture Index (NDMI).These indices have been widely used as comprehensive indicators of the interaction between land cover and solar radiation in the visible and near-infrared regions of the electromagnetic spectrum [26,[39][40][41].
The direct relationship between terrain variables and forest species composition, tree height growth and other stand variables enables these forest variables to be modelled [42,43].First-and second-order terrain parameters were computed from a Digital Elevation Model (DEM) of the study area with a spatial resolution of 15 m [44].Prior to these calculations, a low pass filter with a 5 × 5 moving window was applied to the original DEM in order to reduce the banding effects present in the original file.The first-order terrain parameters selected were elevation, slope, aspect, transformed aspect, profile curvature, plan curvature and curvature, while the second-order terrain parameters were terrain shape index and wetness index.These parameters are potentially related to key features for forest stand development, such as overall climate characteristics, insolation, evapotranspiration, run-off, infiltration, wind exposure and site productivity [42,45].
The sample plots were geolocated in order to extract the average pixel value with an associated buffer of 25 m for each potential predictor.The pixel data were extracted using R statistical software [46] and the "raster" package.Finally, a database was constructed with the mean biomass values for each Forests 2017, 8, 269 8 of 20 plot; the corrected bands of the Landsat-5 TM sensor (6 bands: 1, 2, 3, 4, 5 and 7), the vegetation indexes (7 indexes) and the terrain variables derived from the DEM (9 variables).

Machine Learning Techniques (MLTs)
We compared the performance of three machine learning techniques for estimating the ABG and carbon at stand level: (i) the non-parametric Support Vector Machine (SVM) for Regression technique; (ii) Regression by Discretization based on Random Forest (RD-RF); and (iii) parametric multiple linear regression (MLR).
SVMs have proven useful in remote sensing of forest environments [47,48].The Shevade et al. [49] modification of the Sequential Minimal Optimization (SMO) with a polynomial kernel and a trade-off parameter value of 1 was used in this study for SVM ensemble.On the other hand, RD-RF employed a classifier (random forest, in this case) on a copy of the data in which the response variable is discretized with equal width.The predicted value is the expected value of the mean class value for each discretized interval (based on the predicted probabilities for each interval).In this study, we used the random forest classification algorithm [50].The success of this technique is based on the use of numerous trees developed with different independent variables that are randomly selected from the complete original set of variables.The number of bins for discretization was fixed at 10 and the number of trees fitted was established at 100.Finally, MLR-the technique most commonly used in this kind of study [51]-was used to select the best set of independent variables for the model.The equation was initially built on all descriptors, and descriptors with the smallest standardized regression coefficients were then removed step-wise from the model until no improvement was observed in the estimate of the average prediction error given by the Akaike information criterion [52].WEKA open source software [53] was used to implement all three techniques used.
Several approaches can be used to test the accuracy of supervised learning algorithms.We used the common method of k-fold cross validation.In this process, the data set is divided into k subsets.
In each application, one of the k subsets is used as the test set and the other k-1 subsets form the training set.Error statistics are calculated across all k trials.This provides a good indication of how well the classifier will perform on unseen data.We used k = 10 and computed several standard performance metrics to calculate the goodness-of-fit statistics for each technique.We compared the performance of the MLTs by using the root mean square error (RMSE), the coefficient of determination (R 2 ) and a paired t-test (corrected) based on Student's t-criterion.Finally, we used the selected MLT to generate above-ground biomass and carbon maps for the study area.

Tree-Level Biomass Equations
Table 2 depicts the parameter estimates and goodness-of-fit statistics of cross-validation of the equations systems fitted for each species and groups of species.All the parameter estimates were significant at α = 0.05.Stem wood biomass and AGB estimates were the most accurate, as indicated by the R 2 and RMSE values, whereas foliage and branch biomass estimates were the least accurate.The values of the goodness-statistics calculated using the residuals of cross-validation indicated that the selected equations performed well.The results of the coefficient of determination indicated that the fitted models explained between 57% and 98% of the observed biomass variance per component for all species, with mean values of 91%, 82%, 83% and 76% of the observed variance in stem wood and stem bark, branch and foliage biomass, respectively, whereas the total AGB equations explained on average 93% of the total observed variance in AGB.Nevertheless, as mentioned above, the explained variation was lower than 70% for three of the components: branch biomass in P. engelmannii (59%) and P. leiophylla (63%); foliage biomass in P. herrerae (60%), Q. crassifolia (64%) and Juniperus depeanna (67%); and stem bark biomass in Juniperus depeanna (57%).12 n = number of trees, d = diameter at breast height (cm), h = total height (m), W w = stem wood biomass (kg tree −1 ), W b = stem bark biomass (kg tree −1 ), W br = wood plus bark biomass of branches (kg tree −1 ), W f = foliage (leaves/needles) biomass (kg tree −1 ), W t = total above-ground biomass (W w + W b + W br + W f ) (kg tree −1 ), RMSE, root mean square error.
The root mean square error (RMSE) for the cross-validation ranged between 22.81 and 140.91 kg for total above-ground biomass (W t ) (Table 2), whereas the values of these statistics for stem wood (W w ), stem bark (W b ) branch (W br ) and foliage (W f ) biomass estimates were 15.11-127.6,0.85-27.87,3.02-67.38and 1.72-16.67kg, respectively.The RMSE of W t for all pine and all oak species were 83.04 kg and 134.12 kg, respectively.
The allometric equation (Equation ( 1)) produced the best estimates of all biomass components of P. herrerae, P. lumholtzii, P. leiophylla, P. douglasiana, P. michoacana and P. oocarpa (except stem wood biomass), as well as for the branch and foliage biomass of P. leiophylla, foliage biomass of P. teocote, stem bark biomass of Q. rugosa, and stem bark and foliage biomass of all of the oak species.
A strong relationship between biomass components and total height was found for the groups of pine and oak species (Figure 2).The biomass equation including total height (h) as a predictor was the best equation for all biomass components of P. cooperi, P. durangensis, P. engelmannii, P. teocote (except foliage biomass), P. strobiformis, Arbutus bicolor, Q. sideroxyla, Q. rugosa (except stem bark biomass) and Q. durifolia, as well as for all of the pine species and for stem wood and branch biomass of all of the oak species.AGB and all biomass components of Juniperus depeanna, P. leiophylla and Q. crassifolia were best estimated using the combination of diameter at breast height and total height (d 2 h) (Equation ( 3)). the best equation for all biomass components of P. cooperi, P. durangensis, P. engelmannii, P. teocote (except foliage biomass), P. strobiformis, Arbutus bicolor, Q. sideroxyla, Q. rugosa (except stem bark biomass) and Q. durifolia, as well as for all of the pine species and for stem wood and branch biomass of all of the oak species.AGB and all biomass components of Juniperus depeanna, P. leiophylla and Q. crassifolia were best estimated using the combination of diameter at breast height and total height (d 2 h) (Equation ( 3)).The biomass equation including total height (h) as a predictor was the best equation for all biomass components of P. cooperi, P. durangensis, P. engelmannii, P. teocote (except foliage biomass), P. strobiformis, Arbutus bicolor, Q. sideroxyla, Q. rugosa (except stem bark biomass) and Q. durifolia, as well as for all of the pine species and for stem wood and branch biomass of all of the oak species.AGB and all biomass components of Juniperus depeanna, P. leiophylla and Q. crassifolia were best estimated using the combination of diameter at breast height and total height (d 2 h) (Equation ( 3)).

Above-Ground Biomass Allocation
Biomass distribution varied among tree components and between species.For all species, the largest amount of biomass was located in the stem wood (Figure 3).Stem wood biomass of pine ranged from 63% (P.oocarpa) to 80.9% (P.lumholtzii), while for oak species the values varied between 59.5% (Q.rugosa) and 65.3% (Q.durifolia).General equations for all pine and all oak species estimated stem wood biomass values of 92% and 78%, respectively.
The lowest values for stem bark, branches and foliage biomass were respectively 2.4% (J.depeanna), 11.6% (P.lumholtzii) and 2.4% (P.herrerae), and the highest values for these components were 18.7% (Q.rugosa), 30.5% (A.bicolor) and 7.1% (P.oocarpa), respectively.For all pine species pooled together, the highest quantity of biomass was contained in stem wood (75.4%); followed by branches (14.9%), stem bark (6.4%) and foliage (3.3%).For all oak species, the biomass values for the different components were respectively 62.9%, 22.2%, 10.9% and 4.0% (Figure 3).The biomass equation including total height (h) as a predictor was the best equation for all biomass components of P. cooperi, P. durangensis, P. engelmannii, P. teocote (except foliage biomass), P. strobiformis, Arbutus bicolor, Q. sideroxyla, Q. rugosa (except stem bark biomass) and Q. durifolia, as well as for all of the pine species and for stem wood and branch biomass of all of the oak species.AGB and all biomass components of Juniperus depeanna, P. leiophylla and Q. crassifolia were best estimated using the combination of diameter at breast height and total height (d 2 h) (Equation ( 3)).

Above-Ground Biomass Allocation
Biomass distribution varied among tree components and between species.For all species, the largest amount of biomass was located in the stem wood (Figure 3).Stem wood biomass of pine ranged from 63% (P.oocarpa) to 80.9% (P.lumholtzii), while for oak species the values varied between 59.5% (Q.rugosa) and 65.3% (Q.durifolia).General equations for all pine and all oak species estimated stem wood biomass values of 92% and 78%, respectively.
The proportion of stem biomass in pine trees increased with d (exceeding 50% of AGB for all diameters); although the proportion of stem bark and foliage biomass decreased, it was fairly stable for branch biomass (Figure 4).In all of the oak species, stem wood and branch biomass increased from 0.25 and 0.14 for d = 5 cm to 0.42 and 0.23 for d > 50 cm, respectively.Stem bark of oaks decreased with d from 0.16 to about 0.10, while foliage biomass decreased from 0.13 to 0.03.

Carbon Fractions in Different Tree Components
The mean values of carbon fraction of the pine species were in all cases close to the value provided by IPCC [54] (0.5), except for oak species, in which the average proportion of carbon was 0.45.Sample size was 27 for all species (Table 3).

Carbon Fractions in Different Tree Components
The mean values of carbon fraction of the pine species were in all cases close to the value provided by IPCC [54] (0.5), except for oak species, in which the average proportion of carbon was 0.45.Sample size was 27 for all species (Table 3).

Carbon Fractions in Different Tree Components
The mean values of carbon fraction of the pine species were in all cases close to the value provided by IPCC [54] (0.5), except for oak species, in which the average proportion of carbon was 0.45.Sample size was 27 for all species (Table 3).where Total: total carbon content, s.d = standard deviation.

Biomass and Carbon Estimates in the Permanent Research Plots
The fitted equations were used to estimate the stand level total AGB and C contents in the 429 research plots.The total AGB ranged from 11.06 to 469.42 Mg ha −1 with a mean value of 129.84 Mg ha −1 and C contents ranged from 5.12 to 232.94 Mg ha −1 with a mean value of 63.80 Mg ha −1 .The goodness-of-fit statistics obtained by the different machine learning methods used to model AGB and carbon contents of the 429 sample plots are shown in Table 4.The variables that were most important for estimating AGB and C contents by the MLT evaluated correspond to band 1 and band 7, the vegetation indices NDMI, NBR2 and EVI, and the terrain variables aspect and wetness index.The best AGB estimates were obtained by using Support Vector Machine for Regression (SVM), whereas the best estimates of carbon content were obtained using Regression by Discretization based on Random Forest (RD-RF).However, the results of the paired t-test (corrected) based on Student's t-criterion did not indicate any significant differences in the results produced by the three approaches considered (α = 0.05).A map of AGB and carbon distributions was elaborated using the SVM approach (Figure 5).
The best AGB estimates were obtained by using Support Vector Machine for Regression (SVM), whereas the best estimates of carbon content were obtained using Regression by Discretization based on Random Forest (RD-RF).However, the results of the paired t-test (corrected) based on Student's t-criterion did not indicate any significant differences in the results produced by the three approaches considered (α = 0.05).A map of AGB and carbon distributions was elaborated using the SVM approach (Figure 5).

The Allometric Equations
The equations were based on three model formulations (Equations ( 1)-( 3)) used in previous studies [2,32,[55][56][57].Diameter at breast height was the main explanatory tree variable used to estimate the tree biomass components for all species.Several authors have noted that inclusion of

The Allometric Equations
The equations were based on three model formulations (Equations ( 1)-(3)) used in previous studies [2,32,[55][56][57].Diameter at breast height was the main explanatory tree variable used to estimate the tree biomass components for all species.Several authors have noted that inclusion of total height does not usually lead to a substantial increase in the predictive ability of diameter-based biomass regressions, and they also assume that d is sufficient to yield a reliable tree biomass prediction [58][59][60].However, other authors have found a significant improvement in model fits [61,62] or an increase in the accuracy of the biomass estimates [57,63] when h is also used as a predictor.In this study, d by itself was a good predictor of biomass, but the addition of h as a second variable improved the predictions for several species.A strong relationship between biomass components and total height was found for the groups of pine and oak species, and total height was a significant predictor of all tree biomass components and AGB of 12 out of 17 species and for the group of all pine species, as well as for stem wood and branch biomass of the oak species.The d and h based system of equations yielded poor fits for branch and foliage biomass of P. oocarpa and P. teocote.Lambert et al. [64], Bi et al. [65] and Zhao et al. [66] also found that inclusion of tree height improved the accuracy of predicting the stem biomass but not the crown (needles and branches) biomass components.
Feldpausch et al. [67] reported that the mean relative error in biomass estimates when h was included was 50% lower than when h was excluded.Similarly, we found that inclusion of total height in the biomass equations increased the accuracy of the biomass estimations (measured as RMSE) by between 10.4% (P.teocote) and 53.9% (P.durangensis) for stem wood biomass, by between 2.9% (A.bicolor) and 28.1% (P.strobiformis) for stem bark biomass, by between 13.1% (Q.durifolia) and 20.6% (P.strobiformis) for branch biomass, by between 1.1% (P.engelmannii) and 19.7% (Q.sideroxyla) for foliage biomass and by between 26.7% (Q.rugosa) and 55.9% (P.cooperi) for total tree AGB.These results are consistent with those of António et al. [55], who found that the use of height implied a decrease in the sum of squares of residuals by 72%, 8%, 12% and 10%, for stem wood, stem bark, leaves and branches, respectively.They are also consistent with the findings of Li and Zhao [68], who reported that inclusion of height improved model performance of the fitted equation especially for total above-ground biomass, stem wood biomass and stem bark biomass.
Branch and foliage biomass components are always difficult to estimate with the same accuracy as bole or total AGB.This is illustrated by the fact that inclusion of height in the model only improved the accuracy of branch and foliage biomass estimates by up to 21%.We assume, as Ter-Mikaelian and Parker [69] found, that the influence on height of stand density (i.e., competition), stand structure and site quality may account for some of this variation, and that these variables represent a large proportion of the small contribution of h to explaining the variance in these biomass components for some of the species evaluated.Nevertheless, this was not tested in the present study because the data does not enable discrimination of the effects of these factors or the different silvicultural treatments on the crown structure.
The combined variable d 2 h is usually used in biomass equations [70][71][72], and the accuracy of biomass estimation has been found to increase significantly (measured as R 2 ) when h or d 2 h was also included, in addition to d.In fact, tree biomass is closely correlated with d 2 h, as shown by Parresol [73] and Carvalho and Parresol [74].In the models developed by Parresol [73], height was a good predictor of stem wood but not of stem bark biomass, whereas Carvalho and Parresol [74] obtained the best estimates for stem, crown and total tree biomass of Pyrenean oak by including the variable d 2 h as the sole independent variable in the equation.Bi et al. [65] reached a similar conclusion, reporting that d 2 h performed better for predicting stem and bark components than diameter alone but not for branch and leaf components.In the present study, d 2 h yielded the best estimates of all biomass components and AGB for P. leiophylla, J. depeanna and Q. crassifolia.The results reported here suggest that the best equations for biomass estimation for most species are based on d and h; it is therefore possible to use the biomass equation systems developed for a specific species across different sites in the temperate forests of Durango, assuming that total height is included in the models and that the addition of this variable may take into account different levels of competition induced by different stand density conditions [55].
We found differences in biomass estimates yielded by the equations developed in this study and previously published biomass functions for nine species in northwest Mexico [8].The predictive powers of the new biomass equations are higher than those of previous models developed by Návar [8,9], as noted for the evaluation using the RMSE.For the stem wood biomass component of all pine species pooled together, the RMSE of the new equation was lower (72.12 kg tree −1 ) than the value obtained for the previous equation (105.00 kg tree −1 ); the same behaviour was observed for the oak group, for which the RMSE values were 94.25 and 102.9 kg tree −1 .For the total above-ground biomass, previous equations showed RMSE values of 131.8 and 156.96 kg tree −1 for pine and oaks, respectively; while the equations developed in this study yielded RMSE values of 83.04 and 134.12 kg tree −1 , for these same groups of tree species.Differences relative to the biomass equations proposed by Návar [8,9] may be explained by the following: (i) the previous biomass equations were fitted using tree diameter as the only regressor; however, the inclusion of total height in allometric equations greatly improves the accuracy of individual tree biomass estimations [67,70,75,76], especially in mixed and uneven-aged stands, such as the stands analysed in this study.For example, different species, ages, sizes, crown types and tolerance levels co-exist in the stands under study and biomass equations based only on d cannot therefore describe all possible height-diameter relationships; and (ii) the biomass equations developed by Návar [8] are not additive, and additivity is one of the most important properties of biomass equations [31,74].

Contribution of Components to Total AGB
For all species evaluated, stem wood contributed most to the total AGB, followed by branches, stem bark and foliage.The fact that the highest proportion of the AGB is allocated in the stem has been widely documented, with the proportions ranging from 50% to 92% for different species [4,77,78].Our results are within the same range, varying between 59.5% (Q.rugosa) and 80.9% (P.lumholtzii).
The biomass allocation patterns observed are consistent with those reported by Blujdea et al. [56] for broadleaf species in Romania, Wirt et al. [4] for Norway spruce in Central Europe, Correia et al. [79] for P. pinea L. in Portugal, and by Pajtík et al. [80] for Picea abies [L.] Karst in Slovakia.
As expected, the proportion of biomass allocated to stem wood increased with increasing tree size, whereas the relative contribution of stem bark and foliage to AGB decreased.This is consistent with findings reported by Blujdea et al. [56] and Johansson [58].The pattern was similar for pines and oaks, except for branch biomass.In oaks, branch biomass increased with d, whereas for pines it remained quite constant (around 15%) over the entire diameter range.The decrease in foliage biomass with increasing d can be attributed to competition for light among the trees; small branches and a smaller quantity of leaves occur in larger trees or those trees with a dominant position [81].The foliage biomass of oaks decreased to below 3.5%, while for pines it decreased to about 2%.

Machine Learning Techniques (MLT)
The variables used to estimate AGB with the different machine learning algorithms evaluated can be classified into 2 different groups.In order of decreasing importance, the first group comprises the spectral bands (bands 1 and 7, are the most relevant) and the spectral indices (NDMI, NBR2 and EVI), and the second group comprises the first-and second-order terrain topographical variables (aspect and wetness index) derived from the DEM.Several studies have demonstrated that such variables are usually good predictors of AGB [26,82,83].However, the results may also differ depending on the machine learning method used [26,[84][85][86][87].
The SVM technique yielded the best fits, thus confirming that this type of technique is of great potential for improving biomass estimation, independently of the type of sensor to which it is applied, as demonstrated in recent studies [88].However, no significant differences were found between this method and the other two approaches used.
The values of the goodness of fit statistics for the SVM technique obtained in the present study were slightly lower than those reported by Guo et al. [89], who estimated AGB in Picea crassifolia forests in China using Landsat TM data and two non-parametric methods (k-Nearest Neighbour kNN: R 2 = 0.54 and RMSE = 26.62Mg ha −1 and SVM: R 2 = 0.51 and RMSE = 27.45Mg ha −1 ).They were also lower than those reported by Tian et al. [90] (kNN: R 2 = 0.59) and those obtained by López-Serrano et al. [26] for mixed and uneven-aged forests in the Sierra Madre Occidental of Mexico (R 2 = 0.73 and RMSE = 22.59 Mg ha −1 ).

Conclusions
This study describes a set of simultaneous species-specific allometric equation systems for seventeen temperate species in north-western Mexico.Breast height diameter-only based equations proved best for predicting all biomass components of only four species.Including h or d 2 h as additional predictor in the equation systems only slightly improved stem wood, stem bark and above-ground biomass estimates, but greatly improved the predictions of branch and foliage biomass.
In all species, most of the biomass was allocated in stem wood, followed by branches, stem bark and foliage.The biomass allocation differed between tree components and among species; however, further study is required to clarify which factors affect the allocation patterns.
The proposed species-specific biomass equations can be applied to tree-level data in forest inventories and may also improve the quality of biomass estimates and verify changes in carbon stocks in the temperate uneven-aged multi-species forests in the study area.This research should be viewed as part of an ongoing process, and further sampling is required to provide data for a wider range of species with potentially different growth forms and biomass allocation patterns in order to improve our ability to estimate AGB and carbon for these forests.
The results of this study regarding the use of machine learning techniques to estimate AGB and C content from remote sensors indicate that moderate resolution sensors, such as the Landsat TM5, are potentially useful for monitoring these important variables in real time and at a low cost, because the spectral data are available free of charge.However, the use of new and more accurate biomass equations, such as the systems of equations fitted in this study, combined with high resolution remote sensing data (e.g., Sentinel satellites) to develop predictive machine learning models will enable more realistic and accurate prediction of biomass and carbon stocks at regional level by model-based inference.

Figure 1 .
Figure 1.Location of the study area.Triangles indicate the location of the 429 permanent sampling plots.

Figure 1 .
Figure 1.Location of the study area.Triangles indicate the location of the 429 permanent sampling plots.

Figure 2 .
Figure 2. Relationships between the total height (m) and the tree biomass components stem wood (W w ), stem bark (W b ), branches (W br ) and foliage (W f ) (kg tree −1 ).

ForestsFigure 3 .
Figure 3.Estimated relative contributions of the biomass components to total tree above-ground biomass (AGB) for pine and oak species in Durango state, Mexico.

Figure 4 .
Figure 4. Tree-level biomass allocation patterns of all pine and oak species as a function of individual tree size (d), estimated by the proposed biomass equation systems.

Figure 3 .Figure 3 .
Figure 3.Estimated relative contributions of the biomass components to total tree above-ground biomass (AGB) for pine and oak species in Durango state, Mexico.

Figure 4 .
Figure 4. Tree-level biomass allocation patterns of all pine and oak species as a function of individual tree size (d), estimated by the proposed biomass equation systems.

Figure 4 .
Figure 4. Tree-level biomass allocation patterns of all pine and oak species as a function of individual tree size (d), estimated by the proposed biomass equation systems.

Figure 5 .
Figure 5. Spatial distribution of total above-ground biomass (AGB) (a) and carbon content (b) in the Durango, Mexico forest study area estimated by the Support Vector Machine for Regression (SVM) approach.

Figure 5 .
Figure 5. Spatial distribution of total above-ground biomass (AGB) (a) and carbon content (b) in the Durango, Mexico forest study area estimated by the Support Vector Machine for Regression (SVM) approach.

Table 1 .
Summary statistics of the felled trees and the main above-ground biomass components used to

Table 2 .
Estimated parameters and goodness of fit statistics for the biomass equations for temperate forest species in Durango state, Mexico.

Table 3 .
Specific mean carbon contents of the main tree components of pine and oak species in the study (standard deviation is given in brackets).

Table 3 .
Specific mean carbon contents of the main tree components of pine and oak species in the study (standard deviation is given in brackets).

Table 3 .
Specific mean carbon contents of the main tree components of pine and oak species in the study (standard deviation is given in brackets).

Table 4 .
Summary of the goodness-of-fit statistics for each machine learning method using 10-fold cross-validations for AGB and carbon content.The best results are highlighted in bold type.