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

This paper presents new above-ground biomass (AGB) and biomass components equations for 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. Generalized method of moments was used to simultaneously fit systems of equations for biomass components and AGB, to ensure additivity. Additionally, the carbon content of each tree component was calculated by the dry combustion method, in a TOC analyser. The fitted equations accounted for on average 91, 83, 84 and 78% 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 h or d2h as additional predictor in the d-only based equations systems slightly improved estimates of stem wood, stem bark and total above-ground biomass, and greatly improved the estimates produced by the branch and foliage biomass equations. The fitted equations were used to estimate AGB stocks at stand level from a database on growing stock from 429 permanent sampling plots. Three machine-learning techniques were used to model the estimated stand level AGB and carbon contents; the selected models were applied to map the AGB and carbon distributions in the study area, which yielded mean values of 129.84 Mg ha-1 and 63.80 Mg ha-1, respectively.


Introduction
Better knowledge of carbon stocks and fluxes is needed 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 both practical forestry purposes and 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 obtained in forest inventories [3].In both cases, generic biomass functions are used to quantify the carbon in forests [4] because they improve the accuracy in carbon accounting systems and thus allow accurate planning of whole tree and residual biomass utilization for bioenergy production.
Temperate forests occupy 32 330 508 ha in Mexico, which represents 17 per cent of the country.These are the richest ecosystems in Mexico with some 7000 plant species [5], of which ~150 are species of pines and 170 of oaks, and these represent over 50 per cent of all known pine and oak species [6].Durango is the most important forestry state in Mexico, with 10.5 million ha of forest cover.Timber production in the state is 1.9 million m 3 per year, which represents 32.8% of the national forest production [7].
Mexico is one of the largest emitters of CO2 by deforestation and contributes with 1.6 per cent of global emissions [8], mainly from temperate and tropical forests (12.9 and 54.1 Mt C year −1 , respectively) [9].Thus, most of the states in Mexico are implementing action plans for mitigating the effects of climate change and for accessing economic incentives favoring carbon sequestration in forests.The role of temperate forests in Mexico is important because of their potential to accumulate C and emit large amounts of CO2 into the atmosphere; however, the lack of sets of species-level biomass equations for oak and pine species growing in these forests, in addition to equations that incorporate aspects of forest structure that vary significantly at regional scales is required [10], so that species and site-specific biomass equations must therefore be developed.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 AGB and associated carbon stocks in forests is first required [11].
The objectives of the present study were: (i) to develop species-specific systems of additive equations for predicting total above-ground biomass; and (ii) to model the forest biomass and carbon in the temperate forests in north-western Mexico by using remote sensing Landsat-5 TM imagery, terrain parameters and data from permanent research plots.

Study area
The study was conducted in the temperate uneven-aged and multi-species 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 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 [12].

Biomass data
Biomass estimates for stem wood, stem bark, branches and foliage were obtained from 1336 destructively sampled trees of 17 tree 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.The P. douglasiana and P. michoacana trees were considered together in order to improve the model fit.A general equation for each biomass component for all pines and all oaks trees was also developed.
The fieldwork for aboveground biomass measurement included tree selection, carrying out standing measurements, felling trees, collecting dimensional data, cutting and separating the tree components and weighing fresh components of each biomass component on site.Trees were sampled in 5 cm diameter classes, from 5 cm until the maximum diameter found in the area.The number of trees sampled varied from a minimum of 30 trees 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 following biomass components were considered: stem wood, stem bark, branches (including both wood and bark) and foliage (needles/leaves).For each felled tree, stem diameter outside bark was measured at 0.3, 0.6, 1.3 and thereafter every 2.5 m along the stem to take into account variations in moisture content along the stem.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).Foliage was totally separated from the trunk and weighed on an analytical balance (precision 1 g).A disk 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).On the basis of the ratio of dry biomass to fresh biomass, the biomass of each tree component was calculated 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 analyzer.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    were n= number of trees, d=diameter at breast height (cm), h= total height (m), Ww= wood biomass of stem (kg tree -1 ), Wb= stem bark biomass (kg tree -1 ), Wbr= wood plus bark biomass of branches (kg tree -1 ), Wf= foliage (leaves/needles) biomass (kg tree -1 ), Wt= total above-ground biomass (Ww + Wb + Wbr + Wf) (kg tree -1 ).

Basic models
We used three basic model forms as starting points for model selection (eq.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.
A first regression procedure was used to select the definitive tree variables for each biomass component equation, over the linearized version of the models taking natural logarithms.The significance level for entering and maintaining variables in the model was restricted to 0.001 [2].In this first step, the best model for each tree biomass component and species was chosen.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 [13,14].
The species-specific biomass equation system is formulated as follows: + ε t (5) where Wi represents the tree biomass for the i-th component, Wt is the total tree AGB (i.e. the sum of all the tree biomass components), Xj 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 [14].

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 ® [15].This method produces efficient parameter estimates under potentially heteroscedastic conditions, without specifying the nature of the heteroscedasticity [16], and thus avoids estimating the heteroscedastic error variance.Bi et al. [14] 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.The goodness-of-fit of each biomass fraction model was evaluated using the root mean squared error (RMSE) and the coefficient of determination (R 2 ) (equations 6-7).(7) where Yij and Y 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 above-ground biomass equations developed
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 [17] was used to relate stand biomass and carbon stocks to variables obtained from remote sensors.The total AGB in each stand was calculated by applying the developed 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 on the same dates that the SPIFyS were established (2007 to 2011) and 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 https://espa.cr.usgs.gov).Landsat TM5 bands 1, 2, 3, 4, 5 and 7 were used; band 6 was not used, because of its thermal characteristics [19].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 (NBR); 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 [20][21][22][23].
The direct relationship between terrain variables and forest species composition, tree height growth and other stand variables enables these forest variables to be modelled [24,25].First and second order terrain parameters were thus derived from the 5×5 low pass filtered Digital Elevation Model (DEM) of the study area with a spatial resolution of 15 m [26].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 [24,27].
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 [28] and the "raster" package.Finally, a database was constructed with the mean biomass values for each 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 for Regression (SVM) technique, (ii) Regression by Discretization based on Random Forest (RD-RF), and (iii) parametric multiple linear regression (MLR).
MLR is the technique most commonly used in this kind of study [29]; furthermore, this type of model is easy to understand and is widely used in most scientific disciplines.To select the best set of independent variables, the model 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 [30].On the other hand, Support Vector Machines for Regression (SVM), originating from statistical learning theory, have become a subject of intensive study [31], as they enable the user to deal with highly nonlinear problems [32] such as estimating complex forest structures.These models are developed by a set of vectors that minimize the mean error.SVM are robust in generalization, even when the training data are noisy, and are guaranteed to have a unique global solution that is not trapped in multiple local minima [33].SVM have proven to be useful in remote sensing of forest environments [34,35].The Shevade et al. [36] modification of the Sequential Minimal Optimization (SMO) with a polynomial kernel and a trade-off parameter value of 1 was used for SVM ensemble.Finally, Regression by Discretization based on Random Forest (RD-RF) employs a classifier (random forest, in this case) on a copy of the data in which the property/activity value 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 [37].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.WEKA open source software [38] 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 compute 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 squared 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
The parameter estimates and the goodness-of-fit statistics for the tree biomass component equations for each species are presented in Table 2.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 selected equations fitted the data well and generally explained between 70 and 98% of the observed total and per component biomass variation for all species.Nevertheless, as mentioned above, the explained variation was lower than 70% for two of the components: branch biomass in P. engelmannii (60%) and foliage biomass in P. herrerae (61%).The coefficients of determination were usually highest for stem wood biomass and AGB, and lowest for branch and foliage biomass components.The average variance explained by the biomass equations was as follows: stem wood (91% ± 5.3%), branches (84% ± 9.4%), stem bark (83% ± 9.2%) and foliage (78% ± 7.3%), and the variance explained by the AGB equations was 93% (± 3.6%).The RMSE values ranged between 14.5 kg and 126.7 kg (stem wood biomass), 0.77 kg and 27.7 kg (stem bark biomass), 2.99 kg and 48.1 kg (branches biomass), 1.64 kg and 16.2 kg (foliage biomass) and 4.03 kg and 136.8 kg (AGB).The average RMSE values for AGB of all pine and all oak species groups were 71.2 kg (± 39.5 kg) and 74.1 kg (± 23.2 kg), respectively.
The allometric equation (eq. 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) (eq.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 using the combination of diameter at breast height and total height (d 2 h) (eq.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 2).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 produced stem wood biomass values of 75.5 and 62.9%, 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 of 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%.The biomass allocation is highlighted per species group in Figure 4.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.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 specific carbon content (Mg of carbon per Mg of dry matter) in tree biomass components for the species evaluated is shown in Table 3.The mean values of carbon fraction of the pine species were in all cases close to the value provided by IPCC [39] (0.5), except for oak species, in which the average proportion of carbon was 0.45.

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.
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.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 significance differences in the results (α = 0.05) produced by the three approaches considered.

AGB
The variables that provided most information in the estimation of AGB and carbon content correspond to spectral bands 1 and 7; the vegetation indices NDMI, NBR2 and EVI; and the terrain variables aspect and wetness index (Figure 5).However, the results differ depending on the machine learning method used.The results obtained in terms of the absolute error of estimates for each of the 429 samples based on 10-fold cross validation were used for statistical comparison of the techniques.The comparison is summarized in histograms showing the relative positions reached (rankings) for each technique (Figure 6).In qualitative terms, the SVM technique yielded the best results.

The allometric equations
The equations were mainly based on three model formulations (eq.1-3) used in previous studies [2,14,[40][41][42].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 predictive ability precision of diameter-based biomass regressions, and they also assume that d is sufficient to obtain a reliable tree biomass prediction [43][44][45].However, other authors have found a significant improvement in model fits [46,47] or an increase in the accuracy of the biomass estimates [42,48] 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 it 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 branches biomass of the oak species group.
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. [49], Bi et al. [50] and Zhao et al. [51] 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. [52] reported that the mean relative error in biomass estimates when h was included was 50% less that 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. [40] 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 [53] who reported that the 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 [54] 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 among 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 [55][56][57], and it has found that the accuracy of biomass estimation increased 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 [58] and Carvalho and Parresol [59].In the models developed by Parresol [58], height was a good predictor of stem wood but not of stem bark biomass; whereas Carvalho and Parresol [59] obtained the best estimates for stem, crown and total tree biomass of Pyrenean oak including the variable d 2 h as the sole independent variable in the equation; Bi et al. [50] 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 equations systems developed for a specific species across different sites in the temperate forests of Durango, considering 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 [40].
Differences in biomass estimates between the equations developed in this study and previously published biomass functions for the same species in northwest Mexico [60] were found.The differences with respect to the biomass estimates reported by Návar [60] may be explained by the fact that these equations only use d as an explanatory variable.Finally, we fitted the equations for the four components simultaneously with the total AGB, in order to take into account the correlation between the errors of the models and to guarantee additivity.This restriction has not been taken into account by several authors [60] with the consequent inconsistency in total biomass estimates.

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,61,62].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. [41] for broadleaf trees in Rumania, Wirt et al. [4] for Norway spruce in Central Europe, Correia et al. [63] for P. pinea L. in Portugal, and by Pajtík et al. [64] 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. [41] and Johansson [65].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 [66].
The foliage biomass of oaks decreased to below 3.5%, while for pines it decreased to about 2%.

Machine Learning Techniques (MLT)
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.Several studies have demonstrated that such variables are usually good predictors for estimating AGB [23,67,68].
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 [69].However, no significance differences were found between this method and the other two approaches used.
The values of goodness of fit statistics for the SVM technique obtained in the present study were slightly lower than those reported by Guo et al. [70] 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.62 Mg 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. [71] (kNN: R 2 =0.59) and those obtained by López-Serrano et al. [23] 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 improves stem wood, stem bark and above-ground biomass estimates, but greatly improves 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 developed biomass equations can be applied to tree-level data in forest inventories and may also improve the quality of biomass estimates and verify carbon stocks changes in the temperate uneven-age multi-species forests of the study area.They are also essential tools for accurate estimation of forest residues in the development of bioenergy projects.
The results of this study regarding the use of MLT to estimate AGB and C content from remote sensors indicate that moderate resolution sensors, such as the Landsat TM5, are sufficiently reliable and accurate for monitoring these important variables in real time and at a low cost, because the spectral data are available free of charge.

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

Figure 2 .
Figure 2. Scatter plots showing the relationship between tree biomass components (stem wood, stem bark, branches and foliage) and total height (h).

Figure 3 .
Figure 3.Estimated relative contributions of the biomass components to total tree above-ground biomass (AGB).

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

Figure 5 .
Figure 5. Frequency of appearance (importance) of each attribute (Wrapper feature selection) obtained by each machine learning technique using 10-fold cross validations for AGB (upper) and carbon content (lower).

Figure 6 .
Figure 6.Percentage of sample plots achieved by each technique (ranking) for AGB (left) and carbon content (right) estimates based on 10-fold cross validation.

Table 1 .
Summary statistics of the felled trees and the main aboveground biomass components used to develop the biomass equations.

Table 2 .
Estimated parameters and goodness of fit statistics for the biomass equations for temperate forest species.

Table 3 .
Specific mean carbon contents of the main tree components (standard deviation is given in brackets).
Spatial distribution of total AGB (left) and carbon content (right) in the study area estimated by the Support Vector Machine for Regression (SVM) approach.