Stand-level biomass models for predicting C stock for the main Spanish pine species

Ana Aguirre,Miren del Río,Ricardo Ruiz-Peinado and Sonia Condés

Abstract

Keywords:Martonne aridity index,Dry weight biomass,Carbon stock,National Forest Inventory,Peninsular pine forest,Biomass expansion factor

Background

Forests are fundamental in the global carbon cycle,which plays a key role in the global greenhouse gas balance(Alberdi 2015),and therefore in climate change.As part of the strategy to mitigate climate change,forest carbon sinks were included in the Kyoto Protocol in 1998(Breidenich et al.1998)and subsequent resolutions as the Paris Agreements in 2015.In accordance,countries are requested to estimate forest CO2emissions and removals as one of the mechanisms for mitigating climate change.Based on the international demands,some international institutions request periodic reports on forest indicators which are used in global reports.For example,the State of Europe’s Forest 2015(SoEF 2015)or Global Forest Resources Assessment 2020(FRA 2020)request five-yearly information on accumulated carbon in the biomass of woody species or the accumulated carbon in other sources or sinks.Since the development of these international agreements,numerous countries have made efforts to achieve the main objective of mitigating climatic change.In Spain,for example,the Spanish Ministry for Ecological Transition and Demographic Challenge is developing a data base of the national contribution to the European Monitoring and Evaluation Program(EMEP)emission inventory,which includes Land Use,Land-Use Change and Forestry(LULUCF)sector,with the aim of estimating carbon emissions and removals in each land-use category.Furthermore,annually updated greenhouse gas emission data must be provided for the UNFCCC(United Nations Framework Convention on Climate Change)Greenhouse Gas Inventory Data.

Soil and biomass are the most important forest carbon sinks.The carbon present in soils is physically and chemically protected(Davidson and Janssens 2006),although it is more or less stable depending on the type of disturbances suffered and the environmental conditions(Ruiz-Peinado et al.2013;Achat et al.2015;Bravo-Oviedo et al.2015;James and Harrison 2016).Therefore,the carbon that could be returned to the atmosphere from the ecosystem after a disturbance is mainly contained in the aboveground biomass,which accounts for 70%–90%of total forest biomass(Cairns et al.1997).Carbon stocks and carbon sequestration in tree vegetation are usually estimated thorough biomass evaluation as the amount of carbon in woody species is about 50% of their dry weight biomass(Kollmann 1959;Houghton et al.1996).Although species-specific values can be found in the literature,this percentage is recommended by the Intergovernmental Panel on Climate Change(IPCC)if no specific data is available(Eggleston et al.2006).

There are two main approaches to estimating forest carbon:i)using biogeochemical-mechanisms and ii)the statistical empirical approach(Neumann et al.2016).The second method is more common in forestry since it uses inventory data such as that provided by NFI’s(Tomppo et al.2010)and the data required does not need to be as specific as for the biogeochemicalmechanism approach.Through this approach,biomass and carbon estimates can be obtained using allometric biomass functions and/or biomass expansion factors(BEFs).Biomass functions require variables for individual trees and/or stand variables(Dahlhausen et al.2017),while BEFs convert stand volume estimates to stand dry weight biomass(Castedo-Dorado et al.2012).The BEF method is widely used when little data is available,this being one of the methods recommended in the IPCC guidelines(Penman et al.2003).

BEFs,including their generalization of stand biomass functions depending on stand volume,can be affected by environmental conditions and stand characteristics,such as the species composition(Lehtonen et al.2004;Soares and Tomé2004;Lehtonen et al.2007;Petersson et al.2012;Jagodziński et al.2017).Some authors have also pointed to the dependence of the stand biomass-volume relationship on age or stand development stage(Jalkanen et al.2005;Peichl and Arain 2007;Tobin and Nieuwenhuis 2007;Teobaldelli et al.2009;Jagodziński et al.2017).When age data are not available,as is the case in several NFIs,other variables expressing the development stage can be used as a surrogate of age,such as tree size(Soares and Tomé2004;Kassa et al.2017;Jagodziński et al.2020).In addition,site conditions can influence the relationship between stand biomass and stand volume(Soares and Tomé2004).These conditions can be assessed by means of indicators such as site index or dominant height(Houghton et al.2009;Schepaschenko et al.2018)or directly through certain environmental variables(Briggs and Knapp 1995;Stegen et al.2011).

Most of the information on forests at national level currently comes from the National Forest Inventories(NFIs).Consequently,many countries have adapted their NFIs to fulfil international requirements(Tomppo et al.2010;Alberdi et al.2017).As regards carbon stock,NFIs are widely recognized as being appropriate sources of data for estimating these stocks(Brown 2002;Goodale et al.2002;M?kip??et al.2008),especially at large scales(Fang et al.1998;Guo et al.2010).Although most NFIs are carried out periodically,the frequency does not coincide with the international requirements for data on accumulated carbon and biomass stocks(which may be annual).In the case of the Spanish National Forest Inventory(SNFI),the time between two consecutive surveys is longer than that stated in the international requirements for forest statistics reporting.Hence,the forest indicators from SNFI data should be updated annually in order to fulfill the international requirements.Moreover,the time between two consecutive SNFI is approximately 10 years,although it is carried out a province at a time,so not all the Spanish forest area is measured in the same year.Whereas other countries measure a percentage of their NFI plots each year,distributed systematically throughout the country(allowing annual national estimates to be made,albeit with greater uncertainly),the approach used in Spain is to measure all the plots within a given province,which does not allow for annual data(or indicators)to be extrapolated at national level.As a consequence,indicators must be updated in the same year for all provinces in order to estimate carbon at national level in a given year.A possible approach to updating carbon stocks indicators from SNFI data would be to estimate the stand biomass through tree allometric biomass functions(Neumann et al.2016),although this method would require complex individual tree models to update stand information at tree level(tree growth,tree mortality and ingrowth).Given the strong relationship between stand volume and biomass(Fang et al.1998;Lehtonen et al.2004),estimations of biomass could be also made by updating volume stocks from the SNFI and using BEFs.This option has the advantage that stand volume can often be easily updated through growth models(Shortt and Burkhart 1996)or even by remote sensing(McRoberts and Tomppo 2007).

According to Montero and Serrada(2013),the main pine species(Pinus sylvestris L.,Pinus pinea L.,Pinus halepensis Mill.,Pinus nigra Arn.and Pinus pinaster Ait.)occupy around of 30%of the Spanish forest area as dominant species,which is more than 5 million ha,along with almost half a million ha of pine-pine mixtures.Their distribution across the Iberian Peninsula covers a wide range of climatic conditions(Alía et al.2009),with arid conditions being particularly prominent.Thus,aridity was found to influence the maximum stand density and productivity of these pinewoods(Aguirre et al.2018,2019).Furthermore,pine species were those most used in reforestation programs,so these species play a fundamental role in carbon sequestration.According to the Second and Third National Forest Inventories,the five abovementioned species alone account for a carbon stock of around 250×106Mg C(del Río et al.2017),of which more than half corresponds to two of these forest species(P.sylvestris and P.pinaster).

The main objective of this study was to develop dry weight biomass models for pine forests(monospecific and mixed stands)according to stand volume,exploring whether basic BEFs can be improved by including site conditions and stand development stage.We hypothesized that for a given stand volume the stand dry weight biomass increases as site aridity decreases and that it decreases with the stand development stage.Therefore,the specific objectives were to study the dependence of the models on these factors and to assess the biomass expansion factors when varying these variables for the main pine species studied.The biomass models developed will allow carbon estimates to be updated for a given year when no field data from SNFI surveys are available.

Methods

Data

The data used were from two consecutive completed surveys of the SNFI in the Iberian Peninsula,the Second and Third SNFI(SNFI-2 and SNFI-3),which were carried out from 1986 to 1996,and from 1997 to 2007 respectively,except for the provinces of Navarra,Asturias and Cantabria,where the SNFI-2 surveys were carried out using a different methodology.Data from the SNFI-3 and SNFI-4 were used for these provinces,covering the periods from 1998 to 2000 and from 2008 to 2010,respectively.The initial and final surveys are referred to regardless of the provinces considered.The time elapsed between surveys ranges from 7 to 13 years depending on the province.Data from the final SNFI surveys were used to develop dry weight biomass estimates,while data from the initial surveys,together with volume growth models by Aguirre et al.(2019),were used to evaluate model assessment capability.

The SNFI consists of permanent plots located systematically at the intersections of a 1-km squared grid in forest areas.The plots are composed of four concentric circular subplots,in which all trees with breast-height diameter of at least 7.5,12.5,22.5 and 42.5 cm are measured in the subplots with radii of 5,10,15 and 25 m,respectively.Using the appropriate expansion factor for each subplot,stand variables were calculated per species and for the total plot.For further details of the SNFI,see Alberdi et al.(2010).

The target species were five native pine species in the Spanish Iberian Peninsula:Pinus sylvestris(Ps),Pinus pinea(Pp),Pinus halepensis(Ph),Pinus nigra(Pn)and Pinus pinaster(Pt).Plots located in the peninsular pine forests were used;the criterion for selection being that the density of non-target species should not exceed 5%of the maximum capacity(Aguirre et al.2018).The plots used for each species were those in which the proportion of the species by area was greater than 0.1.Additionally,to allow the application of the results to stands where the volume was updated through growth models,only those plots in which silvicultural fellings affected less than 5%of the total basal area were considered,as this was the criterion used for developing the existing volume growth models(Aguirre et al.2019).

Stem volume was calculated for every tree in the plot according to SNFI volume equations developed for each province,species and stem form(Villanueva 2005).The Martin (1982)criteria were used to obtain volume growth.Dry weight biomass for different tree components was calculated at tree level using equations taken from Ruiz-Peinado et al.(2011),who developed biomass models for all the studied species,using diameter at breast height and total tree height as independent variables.Total tree aboveground dry weight biomass was calculated by adding the weight of stem(stem fraction),thick branches(diameter larger than 7 cm),medium branches(diameter between 2 and 7 cm)and thin branches with needles(diameter smaller than 2 cm).Based on tree data and using the appropriate expansion factors for each SNFI subplot,the stand level volume and dry weight biomass were obtained per species and total plot.

To estimate the aridity conditions for each plot used,the annual precipitation(P,in mm)and the mean annual temperature(Tm,in°C)were obtained from raster maps with a one-kilometer resolution developed by Gonzalo Jiménez(2010).These variables were used to obtain the Martonne aridity index(De Martonne 1926),M,calculated as M=P/(Tm+10),in mm·°C?1.M was chosen as an aridity indicator because of its simplicity and recognized influence on volume growth(Vicente-Serrano et al.2006;Führer et al.2011;Aguirre et al.2019)and maximum stand density(Aguirre et al.2018).Hence,M was expected to have a positive influence on dry weight biomass.

Due to the lack of age information for SNFI plots,the development stage had to be estimated through specific indicators.Tree-size related variables are commonly used as surrogates for stand development stage,one such variable being the mean tree volume(vm),which could be used to correct the lack of age information.The vm was calculated as in Eq.1,where V is the volume of the stand in m3·ha?1,and N is the number of the trees per hectare,both referred to the target species(sp).

A summary of the data used to develop the models is shown in Table 1(note that when a target species wasstudied,other pine species could be included within stands).Figure 1 summarizes the methodology that is described in the following sections.

Table 1 Summary of data used to develop dry weight biomass models.Note that plots where a target species,sp,is studied,other pine species could be present

Biomass estimation models by species

Basic biomass models were developed for each species from SNFIFdata in accordance with the structure used by Lehtonen et al.(2004)(Eq.2)to estimate dry weight biomass(W)from stand volume(V)for the target species.The Basic Model was modified by including the effect of aridity,thus,the Martonne aridity index(M)was added to the Basic Model to obtain the so-called Basic M Model(Eq.3).As regards the model structure,following a preliminary study(not shown)it was decided to include the logarithm of this variable to adapt the Basic Model(Eq.2),modifying the‘a(chǎn)’coefficient according to Eq.3.

where,for plot j in province k,W is the dry weight biomass of the target species in Mg·ha?1,V is the volume of the target species in m3·ha?1,M is the Martonne aridity index,in mm·°C?1;and?is the model error.The coefficient a is the fixed effect,while akis the province random effect to avoid possible correlation between plots belonging to the same province,as the measurements in the different provinces were carried out in different years and by different teams.b and m are other coefficients to be estimated:if coefficient m was not significant for a given species or its inclusion did not improve the Basic Model,M was no longer included in the species model.

To determine how the stand development stage influences the relationships between volume and dry weight biomass for each species,the mean tree volume(vm)was included in the models.This variable also multiplies the coefficient‘(a+ak)’(Eq.4),so that if it was not significant,the final model will be equivalent to the basic one.

where,a,ak,b,c1,p1and m were the coefficients to be estimated and vm is the mean tree volume,all variables referring to the target species.

When fitting the biomass models some bias linked to the stem form was detected.Hence,the next step was to test whether it was possible to correct the model bias by adding the shape of the trees by means of the stand form factor(f)(Eq.5).This variable was also added to multiply the coefficient‘(a+ak)’,thus obtaining the Total Model(Eq.6).

Fig.1 Schematic explanation about how to apply the developed model for future projections.SNFIF is the last Spanish National Forest Inventory available,ΔT is the time elapsed between SNFIF and the projection time T,M is the Martonne aridity index,Origin is the naturalness of the stand(plantation or natural stand),dg is the quadratic mean diameter(cm),Ho is the dominant height(m),RD is the relative stand density,p is the proportion of basal area of the species in the stand,VGE is the volume growth efficiency,IV is the volume increment(m3·ha?1·year?1),N is the number of trees per hectare,V is the volume of the stand(m3·ha?1),vm is the mean tree volume,f is the stand form factor,W is the dry weight biomass,and C is the weight of carbon.The subscript“F”refers to the final SNFI,the last available,while“T”refers at projection time T.The variables with the subscript“sp”refer to the target species,variables without the subscript refer to the stand

where f is the stand form factor;V is the stand volume(m3·ha?1);G is the basal area(m2·ha?1);and H is the mean height of the plot(m),all variables referring to the target species.

where a,ak,b,c1,c2,p1,p2and m were the coefficients to be estimated,f is the form factor of the stand and vm is the mean tree volume,all variables referring to the target species.

The model structure was analysed in a preliminary study where each coefficient in the allometric basic model was parametrized in function of M,vm and f,considering linear and non-linear expansions.The final model structure(Eq.6)was selected because its better goodness of fit in terms of AIC,showing also the lowest residuals.

All models(Eqs.2 to 4 and Eq.6)were fitted using non-linear models with the nlme package(Pinheiro et al.2017)from the R software(Team RC 2014).The coefficients were only included if they were statistically significant(p-value<0.05)and their inclusion improved the model in terms of Akaike Information Criterion(AIC)(Akaike 1974).Furthermore,conditional and marginal R2(Cox and Snell 1989;Magee 1990;Nagelkerke 1991)were calculated as a goodness-of-fit statistic using MuMIn library(Barton 2020).Once selected the model with the lowest AIC,and highest marginal and conditional R2,and to check that the improvement achieved is significant,anova tests were made.

Evaluation of biomass estimation models

In order to evaluate the goodness of fit,an analysis of the four developed models(Eqs.2 to 4 and Eq.6)was performed.The mean errors(Eqs.7 to 9),estimated in Mg·ha?1,as well as mean percentage errors(Eqs.10 to 12)in%were calculated for each model of each species.

where ej=Wj?^Wjand epj=（Wj?^Wj）/Wj;^Wjis the estimated values of dry weight biomass for each plot j,Wjthe corresponding observed values for each plot j,both referring to the target species;and n is the number of plots where the species was present.

Carbon predictions at national level

The models developed(Eqs.4 to 6 and Eq.8)provide estimates of dry weight biomass per species,both in monospecific and mixed stands,which could be transformed to carbon stock,considering the specific data of carbon content in wood given by Ibá?ez et al.(2002)for the five studied pine species(Table 2).

To evaluate the prediction capacity of the fitted models at time T when no field data is available,a simulation from the initial SNFI survey(SNFII)was performed at a national scale,assuming that this was the last available survey.

The first step was to obtain the predicted biomass at time T,where all variables are supposed to be unknown for each species,from the four biomass models developed(Eqs.2 to 4 and Eq.6).To apply these models,it was necessary to obtain the values of all independent variables,updated to year T.This procedure was done as follow:

Table 2 Carbon content of wood for the studied species(Ibá?ez et al.2002)

How to estimate carbon stocks at national level when no data is available

In this section,it is explained how to apply the developed models for predicting the carbon stock at time T required,when no data is available.For this,it is necessary to use some variables of the last Spanish National Forest Inventory available(SNFIF),ΔT years before T.

The first step is to estimate the volume growth efficiency of the target species(VGEsp),which can be estimated using Aguirre et al.(2019)models.These models estimate VGE as function on:

–Origin,makes reference to the naturalness of the stand.It was a dummy variable,with value 1 when the stand was a plantation and 0 when the stand comes from natural regeneration.

–dgsp,is the quadratic mean diameter of the target species.

–Ho,is the dominant height of the stand.

–RD,is the relative stand density(Aguirre et al.2018,Eq.S1),and RDspis only considering the target species.

–psp,is the proportion of the species.

–M,is the Martonne aridity index.

With these variables it is possible to estimate VGEspfor each pine species considered,and using its proportion,also volume growth of each species(IVsp)can be estimated.Note that in monospecific stands IVspis equal to IV total.

Results

Biomass estimation models for each species

Table 3 shows the coefficient estimates together with the standard errors and goodness of fit for the four models developed for dry weight biomass of the five species studied(Eqs.2 to 4 and Eq.6).When the Basic Model(Eq.2)was compared with the Basic M Model(Eq.3)it was observed that aridity(M)was significant in three of the five species and in all three cases it resulted in an improvement in the Basic Model,both in terms of AIC and marginal and conditional R2.The species for which M was not significant in the models were Pt and Pp.Among the species for which M was significant,Ps and Ph showed the greatest increase in conditional and marginal R2,while a slightly negative effect was only detected in the case of Pn(Table 3).

0.9654 0.9659 0.9787 0.9789 0.9454 0.9454 0.9484 0.9750 0.9671 0.9680 0.9756 0.9824 0.9787 0.9790 0.9932 0.9937 0.9828 0.9828 0.9828 0.9829 Table3Coefficientsestimated(a,b,m,c1,p1,c2,p2)andstandarderror(inbrackets)formodelsfromEqs.2to4andEq.6,togetherthestandarddeviationoftherandom M.R2C.R2 0.9609 0.9617 0.9758 0.9761 0.9348 0.9348 0.9358 0.9648 0.9495 0.9495 0.9604 0.9784 0.9762 0.9764 0.9903 0.9911 0.9793 0.9793 0.9793 0.9800 AIC 15,309 15,285 14,414 14,396 4184 4184 4155 3769 13,158 13,104 12,554 11,887 10,751 10,735 9133 9027 9360 9360 9360 9354 StdRnd 0.1441 0.1179 0.0446 0.0192 0.1564 0.1564 0.0808 0.0156 0.0900 0.0685 0.0843 0.0852 0.0970 0.1078 0.0487 0.0178 0.0541 0.0541 0.0541 0.0255 andTotalcorrespondtoEqs.2,3,4and6respectively p2?0.1884(0.0341)?1.0004(0.0330)?0.2534(0.0214)c2?0.5534(0.0143)0.1445(0.0415)variable(StdRnd),AkaikeInformationCriterion(AIC)andmarginalandconditional R2(M.R2andC.R2)p1?0.8141(0.0471)?0.7265(0.0467)?0.0868(0.0130)?0.1967(0.0087)0.5947(0.0613)0.8163(0.1005)?0.5930(0.0312)?0.4745(0.0301)0.0162(0.0066)c1 0.0384(0.0060)0.0536(0.0087)0.4144(0.0236)0.1571(0.0149)0.0766(0.0094)0.1282(0.0173)m 0.0868(0.0178)0.1738(0.0225)0.1980(0.0253)0.1429(0.0233)0.0591(0.0141)0.0547(0.0117)?0.0326(0.0079)0.0275(0.0072)0.0395(0.0076)b 0.7953(0.0040)0.7887(0.0041)0.8482(0.0038)0.8460(0.0038)0.8430(0.0086)0.8430(0.0086)0.8575(0.0087)0.8762(0.0061)0.9466(0.0042)0.9365(0.0043)0.9258(0.0038)0.9132(0.0032)0.8905(0.0039)0.8914(0.0039)0.9422(0.0024)0.9363(0.0024)0.8997(0.0037)0.8997(0.0037)0.8997(0.0037)0.9009(0.0040)a 2.7422(0.0645)2.1193(0.1109)1.0769(0.0612)0.4692(0.0302)2.4602(0.1064)2.4602(0.1064)1.0857(0.0514)0.2988(0.0122)1.2790(0.0299)0.9246(0.0472)1.0488(0.0427)1.6377(0.0570)1.8320(0.0415)2.0679(0.0787)1.0489(0.0330)0.4478(0.0162)1.2275(0.0257)1.2275(0.0257)1.2275(0.0257)0.5757(0.0180)Model Basic BasicM vm Total Basic BasicM vm Total Basic BasicM vm Total Basic BasicM vm Total Basic BasicM vm Total sp Ps Pp Ph Pn Pt sp,arethespeciesanalyzed:PsPinussylvestris,PpPinuspinea,PhPinushalepensis,PnPinusnigra,and PtPinuspinaster.NamesofmodelsBasic,BasicM,vm

The estimates obtained for the coefficients c1and p1in the models that include vm indicate the high importance of this variable for estimating biomass weight.Nevertheless,its influence was less in the case of Pt,as reflected by its low p1value(Fig.2c,Table 3).The coefficients can be significant either as exponents or by multiplying the variables,or in both ways.

The bias observed when fitting the models was corrected by including the stand form factor f.When the Total Model and vm Model were compared,the bias correction was more clearly observed in the Ph model,while for Ps and Pt the inclusion of f only had a slight effect(Table 3).

When the estimation errors were analyzed using the different models(Table 4)it was observed that the bias was always less than 0.2 Mg·ha?1,which in relative terms is equivalent to less than 3%.In general,the models overestimated the biomass weight (negative ME),although for Ph and Pp all the fitted models overestimated the biomass,except the Total Model for Pp.In addition,Pn and Pp were the species for which the greatest reduction in RMSE was observed,comparing the Total Model and Basic Model(greater than 4.5%),while this reduction was the lowest for Pt(around 0.06%).

Fig.2 The selected model(Total Model),showing the dry weight biomass estimations for the target species(W,in Mg·ha?1)according to:a volume of the stand for the target species(V,in m3·ha?1);b Martonne aridity index(M,in mm·°C?1);c mean tree volume(vm,in m3 per tree);and d stand form factor(f).The variable represented in each figure on the x axis,ranges from 1%to 99%of its distribution in the data used,while the rest of the variables remain constant and equal to:V=150 m3·ha?1;M=30 mm·°C?1;f=0.5;and vm=0.5 m3 per tree.Species as in Table 3

Table 4 Model errors calculated through Eqs.7 to 12

Having selected the Total Model as the best model to estimate the dry weight biomass for all species,the influence of each independent variable was analyzed.In Fig.2,the variation of dry weight biomass with each variable was presented,assuming the rest of the variables not represented on the axis remain constant.Figure 2a shows a clear positive relationship between dry weight biomass and stand volume,with Pp being the species producing the highest stand biomass for a given volume,although it was very similar to Ph and Pn.If stand volume(V)is considered constant,it is possible to analyze the variation in W with aridity(Fig.2b),observing that for all species where M was included in the model(Ps,Ph and Pn)the relationship was positive,that is,the higher the M value(less aridity),the higher the W value for a given V.Furthermore,the effect of aridity on this biomass-volume relationship varied according to the species,with Ps being the species for which this influence was the greatest(Fig.2b,Table 3).Analyzing the dry weight biomass variation according to vm(Fig.2c),it was observed that the tendency of the relationship between W and vm was similar for Pp,Pn and Ps,that is,the higher the mean tree volume,the lower the W estimated for a given V.An increase in vm,for a constant V,indicates that the stand is composed of a smaller number of larger trees whereas a decrease in vm indicates that the same stand volume comprising a greater number of smaller trees.Figure 2c shows that the vm effect is more evident when trees are smaller,while the relationship tends to be more constant as the size of trees increases.Note that for Pt and Ph,the vm effect was opposite to that for the other studied species,that is,positive.Figure 2c shows this effect clearly for Ph,despite being the species with the lowest range of vm variation,while for Pt,the influence of vm was only slight,despite being one of the species with the highest range of variation of this variable.As regards the stand form factor(f),in general,W decreased as f approached the unit value(Fig.2d),although in the case of Pt there is a very slight positive effect of f.The influence of f on W was not decisive for Ps and Pn,while it was especially important for Pp and Ph.

Biomass expansion factors

According to the fitted models,the BEF,i.e.stand biomass weight/stand volume,is not constant but rather decreases as the stand volume increases.Figure 3 represents the species BEF variation within the interpercentile 5%–95%range of the species stand volume in monospecific stands for the mean and the extreme values of each of the independent variables in the Total Model.For all species,the estimated BEF values generally varied between 0.5 and 1.5 Mg·m?3,and the lowest estimations were found for Pt,for which the BEF values were almost constant and around to 0.75 Mg·m?3.In contrast,the species for which the highest BEF was obtained was Pp,when f or vm had lower values.BEF estimations for this species could reach values of more than 1.5 Mg·m?3for low stand volume.

Figure 3 shows that the BEF of Pt was always lower than 0.9 and was not influenced by M and hardly affected by vm or f.The BEF values presented little variation in the M range distribution for any of the pine species studied,despite being a statistically significant variable.However,it can be seen in Fig.3 that Ps was the species most affected by aridity.In contrast,the BEF variation for different vm values was evident(Fig.3),being the variable that produced the most change in BEFs for Ps and Pn,although it also affected Pp.Highly variable BEFs values can be observed for Pp and Ph within the f range distribution of the species,while for Ps and Pt this relationship was practically insignificant.If the different species are compared,Pn shows more constant BEF values than the other species,regardless of stand volume.

Carbon predictions at national level

The results confirmed that the Total Model was also that which gave the lowest bias when carbon predictions were update to time T in the pine stands across peninsular Spain(Fig.4).This model allowed carbon estimates with lower errors,both in absolute and relative terms,than the rest of the models,despite all the assumptions described,that is,constant values for both the number of trees per hectare and stand form factor in the elapsed interval considered.

In Fig.4,it can be seen that all models produced overestimations of carbon stocks,except the Total Model,which produced the lowest bias,although it slightly underestimated carbon stock.Figure 4 also shows that the inclusion of the f variable scarcely modified the errors(MAE,RMSE,MAPE and RMSPE),although the bias decreased significantly.When the Total Model was used,the RMSE obtained when making carbon stock predictions for the studied pine species in the Iberian Peninsula was less than 20%,which is slightly higher than 9 Mg·ha?1of C.This Total Model resulted in an important reduction in the bias,reaching around 2%.

Fig.3 Variation of biomass expansion factor(BEF),defined as dry weight biomass(W,in Mg·ha?1)estimated from the Total Model,divided by stand volume(V,in m3·ha?1),for different values of:Martonne aridity index(M,in mm·°C?1);stand form factor(f);and mean tree volume(vm,in m3 per tree).The lines are drawn within the inter-percentile 5%–95%range of stand volume distribution.Solid lines represent the mean value of the variable for each species and dashed and dotted lines represent the 5%percentiles,the mean 95%of the variable distribution for each species

Fig.4 Mean errors for carbon estimates at plot level for the studied pine species throughout peninsular Spain according the four studied models.ME,mean error(in Mg·ha?1 of C);MAE,mean absolute error(in Mg·ha?1 of C);RMSE,Root mean square error(in Mg·ha?1 of C);MPE,mean percentage error(in%);MAPE,mean absolute percentage error(in%);RMSE,Root mean square percentage error(in%)

Discussion

The use of BEFs to estimate biomass at stand level provides an interesting alternative for predicting biomass and carbon stocks in forest systems since stand volume(V)is the only variable required.However,the use of traditional BEFs,mainly as constant values and generally obtained for stands under specific conditions,can result in biased biomass estimates if they are applied under different conditions(Di Cosmo et al.2016).These biases can have a significant impact on estimated carbon in the tree layer when large-scale estimates are made,as is the case of national-scale predictions(Zhou et al.2016).In this study,stand biomass models have been developed that include other easily obtained variables as independent variables,in addition to the stand volume.The fitted models allow us to update the carbon stocks in pine forests across mainland Spain for the five species studied using SNFI data.The strong relationship between stand biomass and stand volume(Fang et al.1998)implies that the Basic Model can provide a good first estimate of biomass.This is confirmed by the results obtained as the Basic Model yields good fit statistics.This suggests that,to a certain extent,the stand volume should absorb the effects of other variables,such as the stand age or stand density,as well as environmental conditions(Fang et al.2001;Guo et al.2010;Tang et al.2016).Therefore,in the development of the different models,the structure of the Basic Model was maintained,expanding its coefficients so that if the specific coefficients corresponding to the effects of M,vm and f were not significant,the Basic Model is returned.However,the models improved for all species with the inclusion of the other variables(Tables 3 and 4),reflecting the fact that stands with the same volume can have different structures leading to different biomass.This is observed in the improvement achieved with the Total Model,both with regard to the goodness of fit of the model and the errors(Tables 3 and 4),indicating less biased and more accurate estimates when the stand characteristics and the aridity conditions(M)are included.

The positive relationship found between the aridity index M and the dry biomass W for a given stand volume supports the findings presented by Aguirre et al.(2019),who reported higher productions in less arid conditions.This positive relationship between M and W suggests greater crown development and higher crown biomass for the same volume in less arid conditions.However,it is important to highlight that the individual tree biomass equations used did not consider this type of within-tree variation in the distribution of biomass with site conditions(Ruiz-Peinado et al.2011).Hence,the observed effect of M must be associated with changes in the stand structure.For example,the variation in vm according to the aridity conditions,that is,the stand V is distributed over more trees of smaller size or fewer larger trees according to the aridity of the site,since the proportion of crown biomass with respect to total biomass varies with tree size(Wirth et al.2004;Menéndez-Miguélez et al.2021).This would entail an interaction between the effect of M and the effect of vm in the models,as reflected in the case of Pn,which varies from negative in the basic model with M to positive for the vm Model and Total Model.However,in general,M is not the most important variable to explain the variation in W(Fig.2b),as can also be observed in the small BEF variation for the studied species in relation with M(Fig.3).

The variable vm,as surrogate of the stand development stage,has a different influence on the models for Ph and Pt than for the rest of the species(Fig.2c).The observed pattern for Ps,Pp and Pn indicates that the relationship between W and V,or the BEF,decreases with vm,i.e.as the stage of stand development increases,as has been observed previously in other studies(Lehtonen et al.2004;Teobaldelli et al.2009).This behavior may be caused by differences in the relationship between the components of the trees.For example,Schepaschenko et al.(2018)observed an important decreasing effect of age on the branch and foliar biomass factors.Similarly,Menéndez-Miguélez et al.(2021)analyzed the patterns of crown biomass proportion with respect to total aboveground biomass of the tree as its size develops for the main forest tree species in Spain.These authors found that in the cases of Ps and Pp,this pattern was decreasing;while for Pn and Pt it was constant(the study did not include Ph).These within-tree biomass distributions would validate the patterns found in the Ps,Pp and Pt models,but not the Pn model.However,Ph presents a totally different BEF behavior with the variation in vm.Analyzing the modular values of the different biomass fractions for this species presented in Montero et al.(2005),it can be observed that the proportion of crown biomass in this species increases slightly with the size of the tree,which could explain the opposite pattern observed in this species.However,this difference could also be due to the equations used to calculate the biomass(Ruiz-Peinado et al.2011),since the maximum normal diameter of the biomass sample used in that study was 44 cm,whereas for the Iberian Peninsula as a whole it was as much as 97 cm(Villanueva 2005).Schepaschenko et al.(2018)also reported that the number of branches in low productive,sparse forest is greater than in high productive,dense forests,which may be a cause for the increasing tendency of W in Ph in relation to vm.

The results indicate an improvement in the models with the inclusion of the stand form factor,although the magnitude of the effect caused by this variable,as well as the improvement in the models,were greater for Pp and Ph than for the rest of the species(Fig.2d,Table 3).To estimate the stand volume,diameter at breast height,total height of the tree and its shape are used,according to species and province available models(Villanueva 2005).However,to estimate stand biomass,the equations applied for the different tree components only depend on the species,the diameter at breast height and the total height of the tree,without considering the shape of the tree(Ruiz-Peinado et al.2011).This difference explains the advisability of considering the stand form factor to avoid biases in the estimates,although it also highlights the need to study the dependence of the biomass equations on the different components of the tree according to their shape.In turn,this shape depends on genetic factors,environmental conditions,and stand structure(Cameron and Watson 1999;Brüchert and Gardiner 2006;Lines et al.2012).

The models obtained underline the importance of considering the environmental conditions and the stand structure(size and shape of trees)when expanding the volume of the stand to biomass.If constant BEF values are used for all kinds of conditions,biomass may be underestimated in younger and less productive stands,while for more mature and/or productive stands it may be overestimated(Fang et al.1998;Goodale et al.2002;Yu et al.2014).These authors also highlight the need to further our understanding of the influence of these factors on the individual tree biomass equations.In this regard,Forrester et al.(2017)found that the intraspecific variation in tree biomass depends on the climatic conditions and on the age and characteristics of the stand,such as basal area or density.The components that mostly depended on these variables were leaf and branch biomass,which suggests that it would be advantageous to have more precise equations for these tree components,which would therefore modify the stand biomass estimates.However,the inclusion of other variables in the tree biomass models in order to improve the accuracy would require a large number of destructive samples from trees under different conditions(site conditions,stand characteristics,age...),which would be difficult to obtain in most cases.

The suitability of SNFI data to develop models has been questioned by several authors(álvarez-González et al.2014;McCullagh et al.2017).One of the main disadvantages is the lack of control about environmental conditions,stand age or history of the stand(Vilàet al.2013;Condés et al.2018;Pretzsch et al.2019).Another shortcoming is the lack of differentiation of pine subspecies in the SNFI,like the two subspecies of Pn,salzmanii and nigra,or those of Pt,atlantica and mesogeensis,which could lead to confusing results such as those obtained for Pt,which was the only species for which the Basic Model improved with the inclusion of both variables together,vm and f.This could suggest that the relationship between volume and shape of trees differs according to the subspecies considered.

Through the models developed(Fig.4),it is possible to provide more precise responses to the international requirements in terms of biomass and carbon stocks.Since the most recent SNFI,it has become possible to update the information at a required time.For this purpose,the least favourable situation was assumed,that is,that the only information available was that obtained from the most recent SNFI.However,the main limitation of the models developed is that they are only valid for a short time period,when the assumptions made can be assumed and when both climatic conditions and stand management do not vary(Peng 2000;Condés and McRoberts 2017).If the elapsed time would be too long for assuming that there is not mortality and that the stand form factor does not vary,the basic model could be applied.Furthermore,to achieve more precise updates,natural deaths and silvicultural fellings must be considered using scenario analysis or by estimating of past fellings(Tomter et al.2016).Besides,a proper validation with independent data was not possible due to lack of such data.When the SNFI-4 is finished for all Spanish provinces,it would be interesting to validate the models developed.

Conclusions

The results reveal the importance of considering both,site conditions and stand development stage when developing stand biomass models.The inclusion of site conditions in the models for Ps,Ph and Pn,indicate that aridity conditions modulate the relationship between the dry weight biomass of a stand(W)and its volume(V),while for Pp and Pt this relationship was not influenced.As hypothesized,it was observed that for a lower aridity,the biomass weight and therefore that of carbon are higher for the same stand volume.

Besides,the results reveal the importance of considering both size and form of trees for estimating dry weight biomass,and therefore to estimate carbon stock.As expected,the relationship between dry weight biomass of the stand and its volume decreases when the stand development stage(vm)increases,except for Ph whose behavior is the opposite,and Pt which is hardly affected by vm.However,the inclusion of this variable reduces the ME,MAE and RMSE for all the studied species,which indicates the importance of its consideration in the dry weight biomass estimation.

Abbreviations

NFI:National Forest Inventory;SNFI:Spanish National Forest Inventory;BEFs:Biomass expansion factors;Ps:Pinus sylvestris;Pp:Pinus pinea;Ph:Pinus halepensis;Pn:Pinus nigra;Pt:Pinus pinaster;M:Martonne aridity index;vm:Mean tree volume;W:Dry weight biomass;f:Stand form factor;C:Carbon weight

In the subscripts

sp:Referred to the target species;T:Any time when no field data is available;I:Initial NFIsurvey;F:Final NFI survey

Authors’contributions

Condés,del Río,and Ruiz-Peinado developed the idea,Aguirre and Condés developed the models,Aguirre programmed the models,and all authors wrote the document.All authors critically participated in internal review rounds,read the final manuscript,and approved it.

Funding

This research received no specific grant from any funding agency in the public,commercial,or not-for-profit sectors.

Availability of data and materials

The raw datasets used and/or analyzed during the current study are available from Ministerio para la Transición Ecológica y el Reto Demográfico of the Government of Spain(https://www.mapa.gob.es/es/desarrollo-rural/temas/politica-forestal/inventario-cartografia/inventario-forestal-nacional/default.aspx).

Declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Author details

1Department of Natural Systems and Resources,School of Forest Engineering and Natural Resources,Universidad Politécnica de Madrid,Madrid,Spain.

2INIA,Forest Research Center,Department of Forest Dynamics and Management,Madrid,Spain.3iuFOR,Sustainable Forest Management Research Institute,University of Valladolid and INIA,Valladolid,Spain.

Received:16 December 2020 Accepted:20 April 2021