New models for estimating the carbon sink capacity of Spanish softwood species

Quantifying the carbon balance in forests is one of the main challenges in forest management. Forest carbon stocks are usually estimated indirectly through biomass equations applied to forest inventories, frequently considering different tree biomass components. The aim of this study is to develop systems of equations for predicting tree biomass components for the main forest softwood species in Spain: Abies alba Mill., A. pinsapo Boiss., Juniperus thurifera L., Pinus canariensis Sweet ex Spreng., P. halepensis Mill., P. nigra Arn., P. pinaster Ait., P. pinea L., P. sylvestris L., P. uncinata Mill. For each species, a system of additive biomass models was fitted using seemingly unrelated regression. Diameter at the breast height and total height were used as independent variables. Diameter appears in all component models, while tree height was included in the stem component model of all species and in some branch component equations. Total height was included in order to improve biomass estimations at different sites. These biomass models were compared to previously available equations in order to test their accuracy and it was found that they yielded better fitting statistics in all cases. Moreover, the models fulfil the additivity property. We also developed root:shoot ratios in order to determine the partitioning into aboveground and belowground biomass. A number of differences were found between species, with a minimum of 0.183 for A. alba and a maximum of 0.385 for P. uncinata. The mean value for the softwood species studied was 0.265. Since the Spanish National Forest Inventory (NFI) records species, tree diameter and height of sample trees, these biomass models and ratios can be used to accurately estimate carbon stocks from NFI data.


Introduction
Southern European forests are characterised by a distinctive set of features. They support high levels of biological diversity (both plant and animal) as a result of the survival of many species in southern European refuges during the glacial periods. Furthermore, they have a harsh, unpredictable climate, difficult socioeconomic conditions and have suffered a long history of over-exploitation accompanied by landscape transformations since ancient times. In these forests, the non-marketable products and services they provide are usually more valuable than their direct yields, especially timber production (Scarascia-Mugnozza et al., 2000). Soil and watershed protection, biodiversity, scenic beauty and, increasingly, recreational use, are the main functions covered by these stands, to which carbon sequestration has recently been added in accordance with international agreements on climatic change mitigation (Kyoto Protocol, UNFCCC, EU Forestry Strategy, Ministerial Conference on the Protection of Forests in Europe).
Hence, in the context of this function as mitigators of the effects of climate change, it is important to estimate the quantity of biomass present in forests, to understand the way in which the biomass accumulates and how it is distributed among the different fractions of the tree. This information will provide a basis for further nutrient studies and facilitate research on the use of biomass in energy production (Schlamadinger and Marland, 1996;Clark et al., 2001). The use of forest inventories as a data source allows us to estimate the quantity of carbon fixed in living vegetation. However, depending on the quality and the amount of information provided by the forest inventory, the accuracy of these estimations will vary. National Forest Inventories have provided the basis for several regional and national-level carbon budgets (Dixon et al., 1994;Goodale et al., 2002).
Indirect approaches such as biomass expansion factors (BEF's) or biomass equations applied to forest inventory data (Brown, 2002) are usually used to quantify carbon sequestration in forests. BEF's convert stem volume or stand volume directly into biomass weight estimates, although they vary depending on growth conditions and stand development, particularly on stand age (Lehtonen et al., 2004(Lehtonen et al., , 2007, stand timber volume (Fang et al., 2001) or tree height (Levy et al., 2004). Therefore, more complex biomass models can provide more accurate estimations than BEF's, hence, are more commonly used to obtain forest biomass estimations (IPCC, 2003). Biomass models are built using destructive, highly costly sampling procedures and relate the dry weight of biomass to dendrometric characteristics; in most cases, the diameter at breast height (d) and/or the total height (h) of the tree (Crow and Laidly, 1980;Pardé, 1980). Softwood species play an important role in the Mediterranean Basin forests due to their widespread distribution and their ecological and socio-economic value. The most important softwood species in Spain include: Abies alba Mill. (silver fir), A. pinsapo Boiss. (pinsapo fir), Juniperus thurifera L. (Spanish juniper), Pinus canariensis Sweet ex Spreng. (Canary Islands pine), P. halepensis Mill. (Aleppo pine), P. nigra Arn. (black pine), P. pinaster Ait. (maritime pine), P. pinea L. (stone pine), P. sylvestris L. (Scots pine) and P. uncinata Mill. (mountain pine). These softwood species occupy more than 9.9 million ha in Spain, of which 6.4 million are pure forests (MARM, 2008).
The information available for estimating forest biomass varies from one species to another. A number of studies have dealt with biomass estimation in P. sylvestris and P. pinaster, although much of this relates to Northern Europe in the case of the former (Marklund, 1988;Lehtonen et al., 2004;Muukkonen, 2007) or to the Atlantic range of the species in the case of Pinus pinaster (Lemoine et al., 1986;Montero et al., 1999;Porte et al., 2002;Balboa-Murias et al., 2006). However, less research has been undertaken with regard to P. nigra (Neirynck et al., 1998;Fattorini et al., 2004), P. pinea (Cabanettes andRapp, 1978;Correia et al., 2010) or P. halepensis (Grunzweig et al., 2007). J. thurifera was studied in Morocco for biomass production using non-destructive methods due to its ecological importance (Montes et al., 2000(Montes et al., , 2002. As far as the other softwood species are concerned (A. alba, A. pinsapo, P. canariensis and P. uncinata), biomass production has not been studied in any depth due to their limited distribution, so scarce information exists in this regard. Recently, Montero et al. (2005) fitted a set of biomass models for the main forest species in Spain (including those mentioned above), which allow us to quantify the biomass and carbon sequestration in forest ecosystems. These allometric biomass equations relate different tree biomass components (stem, different size branches, foliage, total aboveground biomass and root system) to tree diameter, although additivity among the component equations was not considered and each component was independently fitted.
Although allometric equations based on dbh provide one of the easiest and most accurate ways to estimate root biomass from forest measurements (Drexhage and Colin, 2001;Le Goff and Ottorini, 2001), it can be useful to determine the root:shoot partitioning of biomass for the purposes of ecological studies or carbon accounting. These ratios can be applied to individual plants or stands at local, regional or landscape level (Mokany et al., 2006). Moreover, National Greenhouse Gas Inventories under the IPCC, generally employ root:shoot ratios to estimate root biomass and specific values are often unavailable.
In order to improve the existing biomass estimations and to offer more precise information on carbon accumulation in Spanish forests, new biomass equations for conifer species have been developed in this study. The use of methods that guarantee the additivity property among tree biomass components provides consistency between total tree and tree component biomass and also ensures greater statistical efficiency (Parresol, 1999(Parresol, , 2001. Furthermore, the inclusion of tree height in biomass equations as an additional predictor variable, could improve the accuracy of the biomass estimations (Ketterings et al., 2001). The objectives of this study were: i) to determine the extent to which the use of additive methods and the inclusion of tree height as an independent variable improve biomass estimations in the studied species; and ii) to analyse the root:shoot partitioning of biomass for the main softwood species in Spain.

Study area
Individual tree biomass data were collected in representative regions across the natural range of the studied species in Spain. Data for Abies alba were collected in the Pyrenean Mountain Range; for A. pinsapo in the Sierra de Grazalema and Sierra de las Nieves (Southern Spain); for Juniperus thurifera in Guadalajara (Central Spain); Pinus canariensis on the island of Tenerife (Canary Islands); P. halepensis was sampled in the Segura Mountain Range (South East Spain); P. nigra in the Iberian Mountain Range; P. pinaster in the Central Mountain Range (Guadalajara, Central Spain) and the Sierra Morena Mountain Range (Ciudad Real, Southern Spain); P. pinea in the Northern Plateau (Central Spain) and Huelva (South-West Spain); P. sylvestris in the Central Range (Madrid and Segovia, Central Spain) and P. uncinata in the Pyrenean Mountain Range (Fig. 1).

Data
For each species, stands were selected in medium quality sites (medium site index) distributed according to age classes. Average trees and growing conditions were chosen for the destructive sample. Because A. pinsapo is a protected species, the sample trees were not as representative as those of other species and it was not possible to collect belowground samples. Trees were sampled by 5 cm diameter classes, starting at 7.5 cm up to the maximum diameter found in the area. The number of sample trees varied from a minimum of 21 trees in the case of J. thurifera and P. uncinata to a maximum of 305 trees for P. sylvestris. For each sample tree, diameter at breast height (1.30 m) (d), total height (h) and crown height (hc) were measured. The minimum diameters sampled were between 6.2 for P. sylvestris and 10.0 for P. canariensis and P. nigra. The maximum diameters ranged from 41.0 cm for P. uncinata to 77.3 cm for P. nigra (Table 1).
Sampled trees were felled and separated into biomass components in the f ield. The biomass components considered were: stem with bark (commercial volume, up to a top diameter of 7 cm), thick branches (diameter larger than 7 cm), medium branches (diameter between 2 and 7 cm), thin branches (branches with a diameter smaller than 2 cm) and needles (Montero et al., 1999). Estimation of root biomass was only undertaken on a few trees per species and diameter class due to the complexity and cost of the work involved (Table 1). The root component was collected using a backhoe, by digging a trench around the stump and extracting all the roots inside this hole. Using this approach, most of the root system was extracted. Fine roots were not captured by this method.
Each component was weighed in the field (fresh weight) and a representative sample (10 kg) of each of them was taken to the laboratory to be oven-dried at 102°C to constant weight. Thus, it was possible to determine moisture content and estimate dry matter. In those cases where the stem could not be weighed in the field, the diameter was measured at meter intervals up the stem in order to determine the log volumes using Biomass models for the main softwood species in Spain 179  Smalian's formula. Dry weight was calculated by applying the basic wood density for the different species (Gutiérrez Oliva and Plaza Pulgar, 1967). Needles were totally separated and weighed in the case of small trees, while a subsample of thin branches with foliage was taken to estimate foliage mass in larger trees. Maximum and minimum aboveground biomass for the conifer species sampled ranged from 5 kg for P. pinaster to 3,368 kg for P. sylvestris and belowground biomass ranged from 3 kg (P. halepensis and P. uncinata) to 1,193 kg (P. sylvestris).

Biomass equations
Different linear and non-linear equations found in biomass literature (Table 2) were tested for relating the weight of the biomass components to tree variables for each species using diameter at breast height (d) and tree height (h) as independent variables.
In a first step, the best model for each component and species was chosen, based on graphical analysis of residuals and fitting statistics (bias and precision), computing the mean residuals (MRES), root mean square error (RMSE), model efficiency (MEF) (equations 1-3) (Gadow et al., 2001) and the Akaike information criterion (AIC) (Akaike, 1974). The biological behaviour of the model was also evaluated in order to choose the best equations. [1] where y i is the observed value, ŷ i is the estimated value, ȳ i is the mean observed value, n is the number of observations and p is the number of parameters in the model.
Afterwards, the best models selected were simultaneously f itted using joint generalized regression (known as Seemingly Unrelated Regressions, SUR) (Zellner, 1962) to make consistent estimates of the different components with non-linear systems (Parresol, 2001;Bi et al., 2004) (system of equations 4). This method takes into account cross-equation error correlation in order to come up with the sum of total aerial biomass components (additivity property). The employment of this technique guarantees that total aboveground biomass will be the sum of the tree component estimations.
Weighted regression was used to avoid heteroscedasticity, frequently present in biomass data. Each observation was weighted by the inverse of its variance to homogenize the variance of residuals. This weighting factor was estimated through a power function of an independent variable as explained by Parresol (2001) and Balboa-Murias et al. (2006).
The possible presence of multicollinearity was verified through the condition number (Myers, 1990). Model fits were performed using the MODEL procedure in the SAS/ETS software (SAS Institute Inc., 2004).
In order to evaluate the predictive accuracy of the system of equations, they were compared with the equations previously proposed by Montero (2005) for these Mediterranean species whose models were fitted separately for each biomass component using log transformed data and OLS regression. The RMSE and MEF ratios (equations 5 and 6) obtained using both systems Ruiz-Peinado et al. / Forest Systems (2011) 20(1), 176-188  of equations were used in the comparison (Bi et al., 2004).

Results
Best models, independently fitted for each component and specie, were later included in the SUR f itting (Table 3). Tree diameter was the variable that showed the best correlation with biomass weight for all tree components. Tree height was also included in all stem Biomass models for the main softwood species in Spain 181 fraction models, although it did not always appear in the models for the other biomass components. Furthermore, most of the stem fraction models had the standard allometric form: where W s is the weight of stem biomass, d is the dbh, h is the total height, a, b and c are parameters of the model. For all species, the belowground biomass fraction depended exclusively on the diameter variable.
Parameters of the models from SUR fitting and statistics for bias and precision (MRES, RMSE and MEF) are shown in Table 3. Every parameter was significant at the 95% confidence level. Predicted values versus observed values for total aboveground biomass did not show presence of bias in the models fitted (Fig. 2). Due to multicollinearity problem, it was not possible to include an individual model for needles in any species, this latter fraction being included in the thin branch component.
Stem models had good MEF values, being higher than 0.87 for these softwood species. Branch models 182 R. Ruiz-Peinado et al. / Forest Systems (2011) 20(1), 176-188 showed more variability, presenting smaller MEF values for A. alba and P. uncinata. Thick branch models for P. canariensis and P. nigra also yielded high bias values (MRES). Regarding root biomass, the model efficiency was always greater than 0.85 except in the cases of J. thurifera and P. uncinata. The bias of the root models was also higher in these species than in the rest. In some cases, the thick branch component was absent in smaller trees as this fraction only appears when trees reach a certain size. Therefore, the available sample size for this component was smaller than for the rest of components and the equation was only suitable for trees with a larger threshold diameter. For this reason, models for the thick branch component presented a restriction based on a threshold diameter, which varied from 22.5 cm for P. pinea and J. thurifera to 37.5 cm for P. sylvestris. In the cases of A. alba, P. pinaster and P. uncinata, which do not normally have a large number of thick branches, the sample size for this component was so limited that it was included in the medium branch component.
Comparison with previous models developed by Montero et al. (2005) through RMSE and MEF ratios, showed that RMSE ratio values are equal to or less than 1.0 in all cases (the RMSE values for the newer models were lower than those for the previous equations) (Table 4). Large differences in RMSE were revealed when this ratio was small, as in the thick branch fraction in P. pinea or the stem fraction ratio in P. pinaster and P. sylvestris, which all exhibited values below 0.5. MEF ratios presented values equal to or higher than 1.0 (the MEF values for the newer models were equal to or higher than those of the previous ones). The improvement in model efficiency was particularly notable for the thick branch component and the thin branch (with needle) component in P. pinea, as well as for the thick branch fraction in P. halepensis and the medium branch component in P. nigra. The partitioning of tree biomass into stem, crown (branches + foliage) and belowground components is shown in Figure 3, using the models developed for an average tree with a dbh of 35 cm and height calculated from the sample. The stem was the biggest fraction in all cases for these conifer species. The contribution of the stem to total tree biomass varied from 39.9% for P. halepensis or 40.0% for J. thurifera, to 65.7% for A. alba and 61.6% for P. pinaster. Crown fraction was also an important biomass component, with maximum values of 38.7% and 36.6% for P. halepensis and J. thurifera respectively, and minimum values of 11.4% in the case of P. uncinata and 14.9% for P. canariensis. Root biomass accounted for between a third and a fifth of the total tree biomass, with maximum values of 29.3% and 28.8 % for P. uncinata and P. sylvestris respectively and minimum values of of 16.2% for A. alba and 18.6% for P. nigra.
The root:shoot ratios for the species analysed in this study are reported in Table 5. The mean ratios between belowground and aboveground biomass varied from 0.18 (A. alba) to 0.38 (P. uncinata), with a mean value 184 R. Ruiz-Peinado et al. / Forest Systems (2011) 20(1), 176-188  of 0.27 for the softwood species studied. Statistical differences were found between species: A. alba presented the lowest value different from the rest of species; the Pinus species showed similar values between them and J. thurifera and P. uncinata exhibited a large variability and presented similar values.

Discussion
The estimation of forest carbon stocks from forest inventories requires the use of accurate and unbiased biomass models. In this study, new biomass estimation models for the main forest species in Spain have been developed in order to improve on the performance of previous models which did not consider the additivity property (Montero et al., 2005); the latter being a desirable attribute for a system of biomass component equations. The use of the nonlinear seemingly unrelated regressions (NSUR) method to fit the system of equations guarantees this property, giving consistency (Kozak, 1970) and reducing the confidence and prediction intervals of the biomass estimations (Parresol, 1999(Parresol, , 2001. The new models included tree height as an independent variable in some components, resulting in improved model fit statistics. Other authors have also reported improvements in biomass estimations where height and diameter (rather than just diameter) dependent models are used (Lambert et al., 2005;Cienciala et al., 2006). By including tree height, information regarding the competitive environment (stand age, site index, density...) is indirectly considered in the model (Wirth et al., 2004). This fact makes the model more general and permits the use of the equation for different sites (Ketterings et al., 2001). Other independent variables like crown length have been tested in other studies for crown biomass estimation (Carvalho and Parresol, 2003;Antonio et al., 2007). However, in order to assure the applicability of the models, we have not considered this variable because it is not available in the Spanish National Forest Inventory and is not usually measured in forest inventories.
As a consequence of considering simultaneous fitting and the use of tree height as a predictor variable, the additive biomass equations presented in this study represent a considerable improvement on those proposed by Montero et al. (2005) for the studied species. The latter were fitted to the same data, but log transformed data were used and each biomass component was fitted separately, using ordinary least squares regression. The improvement in model efficiency reached around 50% for the thick branch component of P. halepensis and P. pinea as well as for the thin branch (with needles) component of P. pinea (Table 4). For total aboveground biomass, the greatest improvements were found for P. pinea.
In the case of the stem biomass component, all the models fitted were non-linear with an allometric expression including dbh and total height. This expression is very similar to those used for volume estimations, yielding high model efficiency values. Many authors have highlighted the suitability of this combination of variables for stem biomass predictions (Bi et al., 2004;Antonio et al., 2007). For the other tree components (branches of different diameter sizes), models could be either linear or non-linear, with dbh and/or tree height appearing by themselves with different parameters or together in various different combinations. The coefficients related to dbh were in most cases positive numbers, showing that biomass increases with diameter. Conversely, coefficients related to tree height were sometimes negative, particularly for crown fractions, indicating that for the same diameter size, taller trees allocate less biomass to the crown due to the processes involved in competition for light (Lambert et al., 2005).
The ability of the model to predict biomass is lower for branch components than for the stem (Table 3). Branch and foliage biomass are more dependent than bole biomass on tree competition and stand density, hence they present greater variability (Cole and Ewel, 2006;Návar, 2009). The high variability in crown biomass displayed by P. pinea trees is due to frequent pruning for fruit and firewood production.
A number of authors have proposed the use of general equations to estimate aboveground biomass per genus or group of species (Pastor et al., 1983/1984) (Schroeder et al., 1997. However, the variability found in allometric ratios and in biomass component equations among the studied species, suggests that separate equations for each species are essential to accurately estimate biomass per fraction of the tree. This information is required to estimate nutrient stocks, biomass amounts for firewood after treatment or to consider different management options in relation to the nutrient cycle and carbon stocks (Balboa-Murias et al., 2006;Cole and Ewel, 2006;Bravo et al., 2008).
Belowground biomass is not generally considered in biomass studies because of the high-cost and difficulty involved in sampling, even though it makes up a significant part of the total tree biomass. In the case of the species studied, the root weight model depends exclusively on tree diameter, as determined by Drexhage and Colin (2001) or Le Goff and Ottorini (2001).
As regards the relationship between belowground and aboveground biomass (stem, branches and needles), Kurz et al. (1996) found differences in the relationships between softwood and hardwood species. However, Cannell (1982) and Cairns et al. (1997) found no significant differences between species groups and reported general root:shoot ratio figures of 0.26 for conifers and 0.25 for deciduous. The figure for conifers in the present study is in accordance with this general value. Levy et al. (2004) found a mean root:shoot ratio for conifers in Great Britain of 0.359. This mean value differs from our result, although the f igures they present for P. sylvestris (0.301) and P. nigra (0.224) are similar to our figures for these species. Correia et al. (2010) reported a root:shoot ratio for P. pinea of 0.30 for trees at low densities. This value is slightly higher than ours, although our samples were collected at higher densities so there may have been greater competition between trees. Although the root:shoot ratio can vary depending on tree size and stand characteristics, these high observed values highlight the importance of the root fraction in Mediterranean forests.
Although we took samples from a broad range of diameters, only a small number of root samples were taken to evaluate belowground biomass and soil types were not checked as part of this study. Therefore, further research into belowground biomass in Mediterranean forests may be necessary.

Conclusions
The use of seemingly unrelated regression and the large sample size employed to fit the models, result in accurate additive biomass equations for the main Spanish Mediterranean forest species and provide a considerable improvement on the existing equations. Given that the Spanish National Forest Inventory (NFI) identifies the forest species and records tree diameter and height of all sample trees per plot, the models developed could be applied to NFI data, allowing ecosystem-wide, regional and national carbon accounts.
Belowground biomass accounts for a significant proportion of total biomass in Mediterranean ecosystems, which must be considered and quantified in order to obtain complete biomass and carbon estimations.
In Mediterranean forests, where wood production is not the main function, these biomass estimates by component are highly useful to define the best forest management practices to be followed and to identify the role of the forest as a carbon sink.