Predicting biomass of hyperdiverse and structurally complex central Amazonian forests

18 Pages·2016·3.67 MB·English

Biogeosciences,13,1553–1570,2016
www.biogeosciences.net/13/1553/2016/
doi:10.5194/bg-13-1553-2016
©Author(s)2016.CCAttribution3.0License.
Predicting biomass of hyperdiverse and structurally complex central
Amazonian forests – a virtual approach using extensive field data
DanielMagnaboscoMarra1,2,3,NiroHiguchi3,SusanE.Trumbore2,GabrielH.P.M.Ribeiro3,JoaquimdosSantos3,
VilanyM.C.Carneiro3,AdrianoJ.N.Lima3,JeffreyQ.Chambers4,RobinsonI.Negrón-Juárez5,
FredericHolzwarth1,BjörnReu1,6,andChristianWirth1,7,8
1AGSpezielleBotanikundFunktionelleBiodiversität,UniversitätLeipzig,Germany
2BiogeochemicalProcessesDepartment,MaxPlanckInstituteforBiogeochemistry,Jena,Germany
3LaboratóriodeManejoFlorestal,InstitutoNacionaldePesquisasdaAmazônia,Manaus,Brazil
4GeographyDepartment,UniversityofCalifornia,Berkeley,USA
5ClimateSciencesDepartment,LawrenceBerkeleyNationalLaboratory,Berkeley,USA
6EscueladeBiología,UniversidadIndustrialdeSantander,Bucaramanga,Colombia
7GermanCentreforIntegrativeBiodiversityResearch(iDiv)Halle-Jena-Leipzig,Leipzig,Germany
8FunctionalBiogeographyFellowGroup,MaxPlanckInstituteforBiogeochemistry,Jena,Germany
Correspondenceto:D.MagnaboscoMarra([email protected])
Received:4August2015–PublishedinBiogeosciencesDiscuss.:18September2015
Revised:23February2016–Accepted:26February2016–Published:11March2016
Abstract. Old-growth forests are subject to substantial elsandanavailablepantropicalmodel,weobservedsystem-
changes in structure and species composition due to the in- aticbiasesrangingfrom−31%(pantropical)to+39%,with
tensificationofhumanactivities,gradualclimatechangeand root-mean-squareerror(RMSE)valuesofupto130Mgha−1
extreme weather events. Trees store ca. 90% of the total (pantropical). Our first and second best models had both
aboveground biomass (AGB) in tropical forests and precise low mean biases (0.8 and 3.9%, respectively) and RMSE
tree biomass estimation models are crucial for management (9.4 and 18.6Mgha−1) when applied over scenarios. Pre-
and conservation. In the central Amazon, predicting AGB dicting biomass correctly at the landscape level in hyperdi-
at large spatial scales is a challenging task due to the het- verseandstructurallycomplextropicalforests,especiallyal-
erogeneity of successional stages, high tree species diver- lowinggoodperformanceatthemarginsofdataavailability
sity and inherent variations in tree allometry and architec- for model construction/calibration, requires the inclusion of
ture.WeparameterizedgenericAGBestimationmodelsap- predictorsthatexpressinherentvariationsinspeciesarchitec-
plicable across species and a wide range of structural and ture.Themodelofinterestshouldcomprisethefloristiccom-
compositionalvariationrelatedtospeciessortingintoheight position andsize-distribution variabilityof thetarget forest,
layersaswellasfrequentnaturaldisturbances.Weused727 implyingthatevengenericglobalorpantropicalbiomasses-
trees (diameter at breast height ≥5cm) from 101 genera timationmodelscanleadtostrongbiases.Reliablebiomass
and at least 135 species harvested in a contiguous forest assessments for the Amazon basin (i.e., secondary forests)
near Manaus, Brazil. Sampling from this data set we as- still depend on the collection of allometric data at the lo-
sembled six scenarios designed to span existing gradients cal/regional scale and forest inventories including species-
in floristic composition and size distribution in order to se- specific attributes, which are often unavailable or estimated
lect models that best predict AGB at the landscape level impreciselyinmostregions.
across successional gradients. We found that good individ-
ual tree model fits do not necessarily translate into reliable
predictions of AGB at the landscape level. When predict-
ingAGB(drymass)overscenariosusingourdifferentmod-
PublishedbyCopernicusPublicationsonbehalfoftheEuropeanGeosciencesUnion.
1554 D.MagnaboscoMarraetal.:PredictingbiomassofcentralAmazonianforests
1 Introduction In temperate and boreal forests, the size, ontogeny and
site variations have been captured by the development of
Allometriesdescribehowrelationshipsbetweendifferentdi- genericspecies-specificbiomassestimationmodels(Wirthet
mensions (e.g., length, surface area and weight) of organ- al.,2004;Wutzleretal.,2008)basedondatafromhundreds
isms change non-proportionally as they grow (Huxley and of individuals from a single tree species. However, this ap-
Teissier, 1936). The lack of proportionality arises from the proach is prohibitive in the tropics where thousands of tree
factthatorganismschangetheirshapewhiletheygrow(i.e., speciescoexist(Sliketal.,2015;terSteegeetal.,2013).In-
the dimensions differ in their relative growth rates). As one stead, the challenge is to develop generic local or regional
important application, allometric relationships can be used formulations that also generalize across species (Higuchi et
torelatesimpledimensionsoftrees(e.g.,diameteratbreast al.,1998;Limaetal.,2012;Nelsonetal.,1999;Saldarriaga
height,DBH,ortreetotalheight,H)todimensionsmorerel- etal.,1998).Ideally,theycontainpredictorvariablesthat(1)
evantforforestmanagersandbasicecologicalresearch,such jointlycapturealargefractionofthevariationinducedbythe
aswoodvolumeorwholetreebiomass(Brownetal.,1989; underlying morphological and anatomical gradients and (2)
Higuchietal.,1998;Saldarriagaetal.,1998). arestilleasyenoughtoobtainormeasure.
Allometric relationships and biomass estimation models The development and application of such generic models
candiffersubstantiallybetweendifferenttreespecies,espe- poseanumberofchallenges.Findingtheappropriatemodel
cially in species-rich regions with a high variation in tree structureandestimatingthemodelparametersrequiresadata
sizes and architectures such as in the tropical rainforests setwithalargenumberofindividualmeasurementscontain-
(Baninetal.,2012;Nelsonetal.,1999;Poorteretal.,2003). ing the variable of interest (here aboveground biomass, or
This variation reflects differences in growth strategy and AGB) and the predictor variables (i.e., DBH, H, species’
life history, such as tree species occupying different strata SGsandWD).Importantly,thedatasetshouldideallycover
whenmature(e.g.,understory,canopyoremergentspecies), all possible real-world combinations of predictor values in
successionalgroups(SGs)(e.g.,pioneerorlight-demanding ordertoavoiderror-proneextrapolationsandunreliablepre-
species, such as Cecropia spp. and Pourouma spp., in con- dictions. However, in multiple regression models, this pre-
trast to late-successional or shade-tolerant species, such as conditionisrarelymet,notevenbylargedesignmatrices.
Cariniana spp. and Dipteryx spp.) or environmental mi- Theultimatepredictionistypicallyatthelandscapelevel,
crosites (Clark and Clark, 1992; King, 1996; Swaine and whichrequiressummingupindividualpredictionsforseveral
Whitmore,1988). thousands of trees varying in size and species assignment.
Important and highly variable architectural attributes of Thelargerthevariationofpredictorvalueswithinastand,the
tropical tree species include stem shape (e.g., slender to higheristhelikelihoodthatextrapolationerrorsoccur.This
stout form), branch form and branching intensity (e.g., pla- calls for a validation at the landscape level, which requires
giotropic,orthotropicandunbranched),crowncontour(e.g., a plot-based harvest method. For obvious reasons, this has
round,elongatedandirregular),crownposition(e.g.,under- rarelybeenattempted(CarvalhoJr.etal.,1995;Chamberset
story,canopyandemergent),maximumDBHandH (Hallé, al.,2001;Higuchietal.,1998;Limaetal.,2012).
1974; Hallé et al., 1978). In addition, there is large varia- Notable effort has already been made to parameterize
tioningrowthrate(thespeedatwhichacertaintreevolume global/pantropical AGB estimation models (Brown et al.,
isfilled)andconsequentlyinwoodanatomyamongspecies 1989; Chave et al., 2005, 2014). Commonly, these mod-
(Bowman et al., 2013; da Silva et al., 2002; Worbes et al., els are derived using several different data sets, each of
2003). Wood density (WD), which is particularly important which is comprised of relatively few trees and species. Al-
for biomass estimation, varies significantly across regions though few opportunities exist to evaluate theses models at
(Muller-Landau, 2004) and can differ between species by thelandscapelevel,theyareusedworldwideindifferentcon-
morethananorderofmagnitude(Chaveetal.,2006).Given texts, sites and across successional stages. For instance, the
these sources of variation, it is not surprising that different pantropicalmodelfromChaveetal.(2005)(DBH+WDas
allometries were reported when comparing species (Nelson predictors)overestimatedbiomasswhentestedagainsttrees
et al., 1999), successional stages (Ribeiro et al., 2014), on- in Gabon (Ngomanda et al., 2014), Peru (Goodman et al.,
togenies (Sterck and Bongers, 1998) and regions (Lima et 2014), Colombia (Alvarez et al., 2012) and Brazil (Lima et
al., 2012).Unfortunately, transferring suchestimation mod- al., 2012), but it also underestimated the AGB in mixed-
elstoothercontexts–otherspecies,sizeranges,lifestages, species Atlantic Forest stands in Brazil (Nogueira Jr. et al.,
sites or successional stages – typically leads to predictions 2014).
thatdeviatestronglyfromobservations,especiallywhenthe The availability of such generic AGB estimation models
samplingdesigndoesnotallowtheselectionofrelevantdata applicable to many species and contexts is particularly im-
forproperestimationoftheparametersofinterest(Gregoire portant for management, ecological and biogeochemical re-
etal.,2016)orwhenpredictorrangesarelimitedorneglected searchintropicalforestlandscapesthatencompassapartic-
(ClarkandKellner,2012;Sileshi,2014). ularlywidespectrumoffloristicandstructuralvariation.For
example, in the central Amazonian terra firme forests, 1ha
Biogeosciences,13,1553–1570,2016 www.biogeosciences.net/13/1553/2016/
D.MagnaboscoMarraetal.:PredictingbiomassofcentralAmazonianforests 1555
of old-growth forest can hold more than 280 tree species pantropical/globalmodelsaresuitableforlandscapebiomass
(DBH≥10cm) (de Oliveira and Mori, 1999) with a wide assessments–undertheassumptionthattheypredictbiomass
range of architectures and anatomies (Braga, 1979; Muller- satisfactorilyoverallsortsoftropicalforesttypesandsucces-
Landau, 2004; Ribeiro et al.,1999). At the landscape scale, sionalstages.
this region encompasses a mosaic of successional stages
promoted by windthrows (Asner, 2013; Chambers et al.,
2 Materialandmethods
2013; Negrón-Juárez et al., 2010; Nelson et al., 1994). Dis-
turbed areas include a diverse set of species representing
2.1 Studysite
the range from new regrowth to adult survivors, thereby in-
cluding different SGs (pioneers, mid- and late-successional
OurstudysiteislocatedattheEstaçãoExperimentaldeSilvi-
species), tree sizes and with a broader range of architec-
culturaTropical(EEST),a21000haresearchreserve(Fig.1)
turesthanold-growthforests(Chambersetal.,2009;Marra
managed by the Laboratório de Manejo Florestal (LMF) of
et al., 2014; Ribeiro et al., 2014). Once floristic composi-
the Brazilian Institute for Amazon Research (INPA), Man-
tionchangesandstructuralgradientsincreasetothisextent, aus,Amazonas,Brazil(2◦56(cid:48)S,60◦26(cid:48)W).Averagedannual
allometry becomes more complex and reliable landscape- temperature in Manaus was 26.7◦C for the 1910–1983 pe-
level biomass estimates rely on well-designed and well-
riod (Chambers et al., 2004). Averaged annual precipitation
testedgenericbiomassmodels.
ca.50kmeastofourstudysitewas2610mmforthe1980–
Wereporthereanoveldatasetof727treesharvestedina
2000 period (da Silva et al., 2003) with annual peaks of up
contiguousterrafirmeforestnearManaus,Brazil.Thisdata
to3450mm(daSilvaetal.,2002).FromJulytoSeptember
set includes biomass measurements from 101 genera and at
there is a distinct dry season with usually less than 100mm
least 135 tree species that vary in architecture and are from
of rain per month. Topography is undulating with elevation
different SGs (pioneers, mid- and late successional). These
ranging from ca. 50 to 140ma.s.l. Soils on upland plateaus
treesspanawiderangeofDBH(from5to85cm),H (from
andtheupperportionsofslopeshavehighclaycontent(Ox-
3.9 to 34.5m) and WD (from 0.348 to 1.000gcm−3). We
isols), while soils on slope bottoms and valleys have high
used this data set to parameterize generic AGB estimation
sand content (Spodosols) and are subject to seasonal flood-
modelsforcentralAmazonianterrafirmeforestsapplicable
ing(Tellesetal.,2003).Incontrasttofloodplains(i.e.,igapó
across species and a wide range of structural and composi-
andvárzea)associatedwithlargeAmazonianrivers(e.g.,Rio
tional variation (i.e., secondary forests), using various sub-
Negro and Rio Amazonas), valleys associated with streams
setsoftheavailablepredictors;i.e.,size(DBHandH),SGs
andlow-orderriverscanbeaffectedbylocalraineventsand
andWD.
thushaveapolymodalandunpredictableflood-pulsepattern
We next evaluated our models, as well as the pantropical
withmanyshortandsporadicinundationsmainlyduringthe
modelfromChaveetal.(2014)atthelandscapelevelusing
rainyseason(Junketal.,2011).
a virtual approach. We created scenarios of simulated 100
The EEST is mainly covered by a contiguous closed
1ha forest plots by assembling subsets of the 727 known-
canopy old-growth terra firme forest with high tree species
biomasstreesinourdataset.Thesescenariosweredesigned
diversity and dense understory (Braga, 1979; Marra et al.,
tospangradientsin(1)floristiccomposition,byassembling
2014). The terra firme forests are among the predominant
stands with specific proportions of pioneer, mid- and late-
foresttypesintheBrazilianAmazon(Braga,1979;Higuchi
successional species, and (2) size distributions of trees. We
et al., 2004) and ca. 93% of the total plant biomass is
comparedtheknownbiomassoftheseforestassemblagesce- storedintreeswithDBH≥5cm(Limaetal.,2012;daSilva,
narios to predictions based on the generic models, with the 2007).Thetreedensity(DBH≥10cm)intheEESTis593±
goalofansweringthefollowingquestions. 28treesha−1 (mean±99% confidence interval) (Marra et
al.,2014).Treeslargerthan100cminDBHarerare(<1in-
1. Which variance modeling approach and combinations
dividual ha−1) and those with DBH>60cm accounted for
ofpredictorsproducedthebestindividualtreebiomass
only 16.7% of the AGB (Vieira et al., 2004). In the study
estimationmodel?
region,treemortalityratescanbeinfluencedbyvariationsin
2. WhichmodelmostreliablypredictedAGBatlandscape topography(Marraetal.,2014;Toledoetal.,2012).Floristic
level,i.e.,acrosssuccessionalgradients? compositionandspeciesdemographycanalsovarywiththe
verticaldistancefromdrainage(Schiettietal.,2013).
Weexpectedthatthebestmodel,theonereducingbothmean
deviation and error of single and landscape-level biomass 2.2 Allometricdata
prediction, would require species-specific variables as well
as an additional parameter allowing the modeling of het- We used data from 727 trees harvested in this region (dos
eroscedastic variance. Our approach and the independence Santos, 1996; da Silva, 2007), each with measured biomass
of our data set allowed us to evaluate whether it is still im- andpredictor variables. Thisdatasetcomprised 101genera
portant to build local/regional models or whether available andatleast135specieswithDBH≥5.0cm(Table1;alldata
www.biogeosciences.net/13/1553/2016/ Biogeosciences,13,1553–1570,2016
1556 D.MagnaboscoMarraetal.:PredictingbiomassofcentralAmazonianforests
Table1.Summaryofthedatasetappliedinthisstudy.Treeswere
harvested in the Estação Experimental de Silvicultura Tropical,
a contiguous terra firme forest reserve near Manaus, Amazonas,
Brazil.
Variables Old-growth Secondaryforest Secondaryforest
forest (23yearsold) (14yearsold)
NT 131 346 250
SR 82 63 51
DBH 5.0–85.0 5.0–37.2 5.0–33.1
H 5.9–34.5 3.9–27.0 9.0–15.5
WD 0.348–0.940 0.389–1.000 0.395–1.000
AGB 8.3–7509.1 5.4–1690.2 7.5–1562.8
Variables:numberoftrees(NT),speciesrichness(SR),diameteratbreastheight
(DBH)(cm),treetotalheight(H)(m),wooddensity(WD)(gcm−3)and
abovegroundbiomass(AGB)(drymassinkilograms).
Trees were harvested at ground level. For each tree, the
DBH(cm),H (m)andfreshmass(kg)wererecordedinthe
Figure1.StudysiteofterrafirmeforestnearManaus,Amazonas, field by using a diameter tape, a meter tape and a mechani-
Brazil. calmetalscale(300kg×200g),respectively.TheDBHwas
measured before, while H was measured after harvesting.
For trees with buttresses or irregular trunk shape, the di-
are given in Table S1 in the Supplement). The trees were
ameter was measured above these parts. Each tree compo-
harvested through the plot-based harvest method in an old-
nent (i.e., stem, branches and leaves) was weighted sepa-
growth forest and in two contiguous secondary forests (14-
rately. For large trees, stems were cut into smaller sections
year-oldregrowthfollowingslashandburnand23-year-old
before weighing. The mass of sawdust was collected and
regrowthfollowingaclearcut)(Fig.1).Ratherthananindi-
weighted together with its respective stem section. Leaves
vidualselection,ourplot-basedmethodreliesontheharvest-
andreproductivematerial,whenavailable,werecollectedto
ingofalltreesfoundinselectedplots.Thismethodallowsfor
allow species identification accordingly to the APGIII sys-
avalid/faithfulrepresentationoftheDBHdistributionofthe
tem (Stevens, 2012). Botanical samples were incorporated
targetforestsandalandscapevalidationofthefittedmodels
intheEESTcollection.Thewatercontentforeachtreewas
(Higuchietal.,1998;Limaetal.,2012).
determinedfromthreediscs(2–5cminthickness),collected
Before selecting plots, we surveyed both the old-growth
fromthetop,middleandbottomofthebole,andsamples(ca.
andsecondaryforeststoassurethatnostrongdifferencesin
2kg) of small branches and leaves. The samples were oven
structure and floristic composition existed and that the se- driedat65◦Ctoconstantdrymass.Thedrymassdatawere
lected patches were representative of our different succes-
calculatedbyusingthecorrespondingwatercontentofeach
sional stages. In the old-growth forest the trees were har-
component(Limaetal.,2012;daSilva,2007).Drymassfor
vested in eight plateau and three valley plots (10m×10m)
each tree was used for subsequent model fits and compar-
randomlyselectedwithinanareaof3.6ha(daSilva,2007).
isons.
In each of the secondary forests the trees were harvested in
five plots (20m×20m), each randomly selected within a
2.3 Species’architectureattributes
1ha plateauarea (dos Santos, 1996; daSilva, 2007). By in-
cludingtreesfromsecondaryforestswewereabletoincrease
Each of our tree species or genera was assigned to one of
the variation in floristic composition and consequently the
three SGs known to vary in their architecture, namely pi-
rangeofspecies-relatedvariationinarchitectureandallome-
oneer, mid- and late-successional groups. To make this as-
try(Table1andTableS1).Sinceoursecondaryforestswere
signment,weconsideredseveralattributesrelatedtospecies’
inserted in the contiguous matrix from which old-growth
architecture (i.e., shape and life history), growth position
plotsweresampled,wealsocontrolledfortheeffectsofim-
(i.e., stratum), morphology, wood density and ecology (Ta-
portant drivers of tree allometry and architecture, such as
bleS1andTableS2).Wevalidatedthisapproachbycheck-
variations in environmental conditions (e.g., soil, precipita-
ing our assignments against those of classic studies (Clark
tionratesanddistribution),foreststructureandwooddensity
and Clark, 1992; Denslow, 1980; Saldarriaga et al., 1998;
(Banin et al., 2012); the last is intrinsically related to varia-
Shugart and West, 1980; Swaine and Whitmore, 1988), lo-
tioninfloristiccomposition.
cal/regionalstudiesconductedintheAmazon(Amaraletal.,
2009; Chambers et al., 2009; Kammesheidt, 2000; Marra
Biogeosciences,13,1553–1570,2016 www.biogeosciences.net/13/1553/2016/
D.MagnaboscoMarraetal.:PredictingbiomassofcentralAmazonianforests 1557
et al., 2014) and species descriptions available in the Mis- dicates significant collinearity among variables (Graham,
souri Botanical Garden (http://www.tropicos.org), species- 2003;Petraitisetal.,1996).Modelseries1–4hadVIF<1.5
Link(http://www.splink.cria.org.br)andListadeEspéciesda (Table2),whichindicatednosignificantcollinearityamong
FloradoBrasil(http://www.floradobrasil.jbrj.gov.br/).More predictors. For model series 5–8, we found VIF>2.0 for
importantly, we considered empirical field observations, ar- DBHandH,whichindicatessignificantcollinearitybetween
chitecturalinformationfromourdatasetanddataforspecies these two variables. This pattern was previously reported
presence/absence from a network of permanent plots repre- for other data sets from Amazon and other tropical regions
sentingawiderangeofsuccessionalstagesincentralAma- (Limaetal.,2012;Ribeiroetal.,2014;Sileshi,2014).
zon (Table S2). This network includes plots in old-growth We fit models representing the eight different predictor
forests(LMFunpublisheddata(censusfrom1996to2012); combinations to our entire data set of 727 trees using three
da Silva et al., 2002), secondary forests (Carvalho Jr. et al., variancemodelingapproaches:nonlinearleastsquare(NLS),
1995;dosSantos,1996),andsmallandlargecanopygaps(≥ ordinaryleastsquarewithlog-linearregression(OLS)anda
ca.2000m2)createdbywindthrowsthatare4,7,14,17,24 nonlinearapproachinwhichwemodeledtheheteroscedastic
and27yearsold(LMFunpublisheddata;Marraetal.,2014). variance of the data set (MOV). In the MOV approach we
Since reported WD values for the same species or gen- modeledthevarianceasafunctionofDBHwithanormally
era can vary strongly among different studies (Chave et al., distributedresidualerror:
2006) and sites (Muller-Landau, 2004), we compiled WD
valuesmainlyfromstudiescarriedoutintheBrazilianAma- εi =N(cid:0)yˆi,σi(cid:1), (1)
zon (Chave et al., 2009; Fearnside, 1997; Laurance et al.,
whereiisthesubscriptforindividuals(i=1,...,n)andσ is
2006; Nogueira et al., 2005, 2007). For species where WD i
modeledwithaheteroscedasticvarianceaccordingto
data were not available for the Brazilian Amazon, we con-
sideredstudiesfromotherAmazonianregions(Chaveetal., σ = c ·DBHc2. (2)
i i i
2009). For species where no published WD was available,
orwheretheidentificationwascarriedouttothegenuslevel Modelseries1(M11,M12andM13)usedDBHasthesole
(63 in total), we used the mean value for all species from predictor(Table2).Formodelseries2(M21,M22andM23),
thesamegenusoccurringincentralAmazon.Fortreesiden- weallowedtheb regressionparametersandc heteroscedas-
tified only to the family level (seven in total), we used the tic variance to vary according to the SG assignment (1, 2
mean value of genera from that family excluding those not or 3). This approach allowed us to account for differences
reportedinthestudyregion(TableS1). amongtheSGswithoutsplittingthedatasetintothreediffer-
entgroups.Thismethodhasincreasedanalyticalpowerand
2.4 Statisticalanalyses allowedustoassesstherelationshipsbetweentreeallometry
andarchitecture.
2.4.1 Individualtreebiomassestimationmodelfits Formodelseries3(M31,M32andM33),weignoredthe
SGassignmentbutintroducedWD(whichdidnotcorrelate
The AGB estimation models we applied varied in the num- stronglywithSG).Formodelseries4(M41,M42andM43)
berandcombinationofourpredictorvariables(eightcombi- weallowedeachSGtohaveitsownwooddensityeffect.For
nations/series) as well as the strategy of modeling the vari- modelseries5and6,wereplacedtheWDwithH.Inmodel
ance (three model types – see below), yielding a set of series5(M51,M52andM53),werestrictedtheSGvariation
24 candidate models (Table 2). We used DBH (cm), WD ofbandc,whileinseries6(M61,M62andM63)weallowed
(gcm−3)andH (m)aspredictors.Furthermore,weusedthe these parameters to vary according to SG. For model series
species’SGassignmentasa“categoricalpredictor”(factor1 7 (M71, M72 and M73), we combined DBH, H and WD
is pioneer, 2 is mid-successional and 3 is late-successional butrestrictedtheSGvariationofbandc.Finally,formodel
species), thereby representing functional diversity along a series 8 (M81, M82 and M83), we combined DBH, H and
main axis of tree successional strategies, functional and ar- WDandallowedbandctovarywithSG(Table2).
chitecturalvariation.Dependingonthemodel-typeparame- In contrast to prior approaches, we did not test mod-
ters, the continuous variables were allowed to vary for cap- els based on compound (e.g., log[AGB] ∼ log[b1] +
turing the successional aspects of functional diversity. We b2[logDBH2HWD]) or quadratic/cubic derivatives (e.g.,
considertheSGgroupingfactorasintegralpartofthemodel. log[AGB] ∼ log[b1] + b2[logDBH] + b3[logDBH2] +
FittingallSGsinonemodelinanMarkovchainMonteCarlo b4[logDBH3])(Brownetal.,1989;Chaveetal.,2005,2014;
context is different than fitting separate models because the Ngomandaetal.,2014).Thesestructureswouldhavelimited
joint model also absorbs the covariance structure of the pa- our ability to include biological variation by defining SG-
rameters across groups, especially in models where not all specificparametersforDBH,H andWDseparately.
parametersareallowedtovaryamongSGs. We fit the AGB estimation models with non-informative
Wetestedvariablesforcollinearitybycalculatingthevari- uniform priors using WinBUGS 1.4.3 (Lunn et al., 2000;
ance inflation factor (VIF). A conservative VIF>2.0 in- Spiegelhalteretal.,2002).Foreachmodel,threechainswere
www.biogeosciences.net/13/1553/2016/ Biogeosciences,13,1553–1570,2016
1558 D.MagnaboscoMarraetal.:PredictingbiomassofcentralAmazonianforests
Table2.Testedequationsforestimatingabovegroundtreebiomass(AGB)inaterrafirmeforestnearManaus,Amazonas,Brazil.
Series Model Equation Variance VIF
modeling(range)
approach
1 M11 AGB∼b1DBHb2 NLS 1
M12 log(AGB)∼log(b1)+b2(logDBH) OLS 1
M13 AGB∼b1DBHb2 MOV
2 M21 AGB∼b1[SG]DBHb2[SG] NLS 1.001
M22 log(AGB)∼log(b1[SG])+b2(logDBH[SG]) OLS 1.005
M23 AGB∼b1[SG]DBHb2[SG] MOV
3 M31 AGB∼b1DBHb2WDb3 NLS 1.007
M32 log(AGB)∼log(b1)+b2(logDBH)+b3(logWD) OLS 1.017
M33 AGB∼b1DBHb2WDb3 MOV
4 M41 AGB∼b1[SG]DBHb2[SG]WDb3[SG] NLS 1.016–1.468
M42 log(AGB)∼log(b1[SG])+b2(logDBH[SG])+b3(logWD[SG]) OLS 1.017–1.395
M43 AGB∼b1[SG]DBHb2[SG]WDb3[SG] MOV
5 M51 AGB∼b1DBHb2Hb3 NLS 3.382
M52 log(AGB)∼log(b1)+b2(logDBH)+b3(logH) OLS 3.342
M53 AGB∼b1DBHb2Hb3 MOV
6 M61 AGB∼b1[SG]DBHb2[SG]Hb3[SG] NLS 1.019–3.439
M62 log(AGB)∼log(b1[SG])+b2(logDBH[SG])+b3(logH[SG]) OLS 1.010–3.360
M63 AGB∼b1[SG]DBHb2[SG]Hb3[SG] MOV
7 M71 AGB∼b1DBHb2Hb3WDb4 NLS 1.014–3.428
M72 log(AGB)∼log(b1)+b2(logDBH)+b3(logH)+b4(logWD) OLS 1.038–3.469
M73 AGB∼b1DBHb2Hb3WDb4 MOV
8 M81 AGB∼b1[SG]DBHb2[SG]Hb3[SG]WDb4[SG] NLS 1.523–3.624
M82 log(AGB)∼log(b1[SG])+b2(logDBH[SG])+b3(logH[SG])+b4(logWD[SG]) OLS 1.422–3.547
M83 AGB∼b1[SG]DBHb2[SG]Hb3[SG]WDb4[SG] MOV
Predictors:diameteratbreastheight(DBH)(cm),species’successionalgroup(SG)(pioneers,mid-andlatesuccessional),treetotalheight(H)(m)andwooddensity(WD)
(gcm−3).Variancemodelingapproach:nonlinearleastsquare(NLS),ordinaryleastsquarewithlog-linearregression(OLS)andnonlinearwithmodeledvariance(MOV).
SinceNLSandMOVrelyonthesameequation,theyhaveanaloguevariationinflationfactorvalues(VIF).
runinparallel,andconvergenceoftheposteriordistribution determination(R2adj)andtherelativestandarderror(S ).
yx%
foreachparameterwasassessedbyconvergenceoftheratio TheS wascalculatedasfollows:
yx%
ofpooledtomeanwithin-chaincentral80%intervalsto1or (cid:18) (cid:19)
2s
bythestabilityofbothintervals(BrooksandGelman,1998; S = √ , (3)
yx%
yˆ N
BrooksandRoberts,1998).
To select the best model we calculated the deviance wheres,yˆandN arethestandarddeviationoftheregression,
information criterion (DIC). The DIC is a generalization the mean of the focal independent variable and the number
of Akaike’s information criterion and consists of a cross- of observations, respectively. As in all allometric data sets
validatorytermthatexpressesboththegoodnessofthefitand relating linear to volume-proportional data, there is indeed
the models’ complexity. The lower the value the higher the heteroscedasticity in our data, which restricts the use of the
predictiveabilityandparsimony(Spiegelhalteretal.,2002). S for model selection. Nonetheless, this measure is pre-
yx%
We also checked whether the 95% credible intervals of the scribedforassessingmodels’uncertainty(IPCC,2006)and
parameter’sposteriordistributionsexcluded0.However,we is commonly used (Chave et al., 2014; Lima et al., 2012;
did not attempt to test the null hypothesis that a particular Ribeiroetal.,2014;Sileshi,2014).
parameteris0(Bolkeretal.,2013;Bolker,2008).Contrasts For the OLS approach including log-transformed vari-
were evaluated by monitoring differences between parame- ables, we calculated the S using untransformed data. To
yx%
ters or predictions based on their posterior distribution. For correctforthebiasintroducedbythelog-transformeddata,a
communicating the results we consider two parameters sig- correctionfactor(CF)wascalculatedasfollows:
nificantlydifferentifthe95%credibleintervaloftheposte- !
SSE2
riordistributionoftheirdifferencedoesnotinclude0. CF=exp , (4)
Toallowforcomparisonsofdifferentmodelstructuresand 2
approaches with the available literature, we calculated the
where SSE is the standard error of the estimate (Sprugel,
coefficientofdetermination(R2),theadjustedcoefficientof
1983).
Biogeosciences,13,1553–1570,2016 www.biogeosciences.net/13/1553/2016/
D.MagnaboscoMarraetal.:PredictingbiomassofcentralAmazonianforests 1559
strainedourfloristiccompositionscenariostoastemdensity
of 1255treesha−1 (DBH≥5cm) typical for the old-growth
terra firme forests at the EEST (Lima et al., 2007; Marra et
al.,2014;Suwaetal.,2012).
To address variations in size distribution, we varied the
proportion of small and big trees fixing a threshold value
of 21cm, which represents the mean DBH (trees with
DBH≥10cm)ofourstudiedforest(Marraetal.,2014).Our
size-distributionscenariosincludedasmall-sizedstand,with
90% from small (DBH<21cm) and 10% from big trees
(DBH≥21cm); a mid-sized stand with equal numbers of
trees smaller and greater than or equal to 21cm in DBH;
anda large-sizedstand,with10%small and90% bigtrees
(Fig.2bandd).Asforourfloristiccompositionscenarios,in
order to produce reliable size-distributions, we constrained
our sampling effort to a basal area value of 30.3m2ha−1
also typical of our studied old-growth forest (trees with
DBH≥5cm) (Marra et al., 2014; Suwa et al., 2012). Both
our floristic and size-distribution scenarios produced the J-
Figure 2. Sampling schemes applied to assemble the six forest
inversedistributionpattern,typicaloftropicalforests(Clark
scenarios designed to reflect changes in floristic composition and
andClark,1992;Denslow,1980).
sizedistributionoftrees,typicalofcentralAmazonianterrafirme
AGBatthelandscapelevelwasdeterminedbyaddingup
forests.
themeasuredAGBfor“sampled”treesineachscenario.To
test how well our biomass estimation models predicted the
2.5 Landscape-levelbiomasspredictionsacross AGB at the stand level, we related biases and root-mean-
scenarios squareerror(RMSE).Inordertoassesstheaccuracyofdif-
ferent predictions in the context of models’ uncertainty, we
To evaluate the models outlined in Table 2, we predicted additionally reported the overall performance of the tested
AGBatthelandscapelevelforsixforestscenariosassembled modelsalongallforestscenarios.Whendoingso,wepresent
byastratifiedrandomselectionofindividualtreesfromour the bias and RMSE in the same unit (Mg), which allow for
data set of 727 trees. Our scenarios were designed to span assessingthemagnitudesofdeviationsinmodelpredictions
a successional gradient created by natural disturbances in (Gregoire et al., 2016; McRoberts and Westfall, 2014). Be-
which the interaction of tree mortality intensity and species cause data on tree height are normally unavailable or esti-
vulnerability and resilience produce complex communities mated imprecisely in Amazonian forest inventories, we fo-
varyinginspeciescompositionandsizedistributionoftrees cused on models including only DBH, WD and SG as pre-
(Chambers et al., 2009, 2013; Marra et al., 2014). We as- dictors(modelseries 1–4).Inadditionto the“internaleval-
sembledthreescenariostoreflectvariationsinfloristiccom- uation”ofourmodels,wetestedthepantropicalmodelfrom
position and three scenarios to reflect variations in size dis- Chaveetal.(2014):
tribution. Each scenario was sampled 100 times, resulting
in100 1haplots perscenario withdifferentcombination of lo(cid:92)gAGB∼−1.803−0.976E+0.976(cid:2)logWD(cid:3)
trees randomly (with replacement) assembled according to
thescenario-specificdesignprinciples. +2.673(cid:2)logDBH(cid:3)−0.0299[logDBH]2, (5)
Toaddresstheeffectofvariationsinfloristiccomposition
on estimated AGB, we created scenarios where we varied which was parameterized with data from 4004 trees
theproportionofpioneer,mid-andlate-successionalspecies. (DBH≥5cm)harvestedin53old-growthandfivesecondary
The early-successional scenario comprised 50% from trees forests.ThismodelhasDBH,H (estimatedfromaDBH:H
sampled randomly from the species classified as pioneer, relationship),WDandavariableE(environmentalstress)as
40% from mid- and 10% from late-successional species predictorsandwassuggestedforestimatingtreeAGBinthe
(as survivors of disturbances). The mid-successional sce- absenceofheightmeasurements.
nario comprised 10% from trees sampled randomly from WeperformedallanalysisusingtheR3.2.1softwareplat-
the species classified as pioneer, 70% from mid- and 20% form(RCoreTeam,2014).WeusetheR2WinBUGS(Sturtz
from late-successional species. The late-successional sce- et al., 2005) package for running WinBUGS from R and
nario comprised 10% from trees sampled randomly from theggplot2package(Wickham,2009)forproducingfigures,
the species classified as pioneer, 40% from mid- and 50% withtheexceptionofFig.1,whichwasproducedintheEnvi-
from late-successional species (Fig. 2a and c). We con- ronmentforVisualizingImagessoftware(ENVI,ITTIndus-
www.biogeosciences.net/13/1553/2016/ Biogeosciences,13,1553–1570,2016
1560 D.MagnaboscoMarraetal.:PredictingbiomassofcentralAmazonianforests
Table3.Statisticsofabovegroundbiomass(AGB)estimationmodelsfitinaterrafirmeforestnearManaus,Amazonas,Brazil.SeeTable2
forpredictorsandappliedvariancemodelingapproachesandTableA3forthemodels’parameters.
Series Model Dev pD DIC R2 R2adj Syx% CF
1 M11 9694.5 2.919 9697.4 0.894 0.894 3.130
M12 6808.0 2.990 6811.3 0.865 0.865 3.542 1.066
M13 6821.0 3.856 6825.2 0.864 0.864 3.544
2 M21 9216.0 3.773 9219.9 0.946 0.945 2.259
M22 6751.0 6.943 6758.3 0.557 0.540 6.458 1.061
M23 6739.0 10.465 6749.5 0.558 0.554 6.381
3 M31 9291.0 4.052 9294.7 0.949 0.939 2.373
M32 6683.0 4.062 6687.0 0.885 0.884 3.280 1.056
M33 6698.0 4.918 6702.5 0.865 0.865 3.527
4 M41 9057.0 2.303 9059.8 0.957 0.956 2.030
M42 6657.0 10.006 6667.5 0.701 0.699 5.215 1.054
M43 6649.0 13.059 6661.6 0.701 0.699 5.239
5 M51 9479.0 0.023 9479.3 0.921 0.921 2.702
M52 6680.0 4.017 6684.3 0.899 0.898 3.060 1.055
M53 6720.0 4.674 6724.7 0.897 0.896 3.103
6 M61 9183.9 −71.746 9112.2 0.948 0.947 2.214
M62 6614.0 10.078 6624.1 0.754 0.750 4.845 1.050
M63 6631.0 11.754 6642.9 0.740 0.737 4.896
7 M71 8998.0 0.951 8999.1 0.959 0.959 1.942
M72 6570.0 5.023 6574.9 0.934 0.933 2.480 1.047
M73 6610.0 5.697 6615.4 0.922 0.920 2.707
8 M81 8812.0 −42.073 8770.3 0.968 0.967 1.719
M82 6548.0 13.031 6561.3 0.811 0.804 4.200 1.046
M83 6566.0 13.778 6580.0 0.801 0.800 4.262
Parameters:models’deviance(Dev),effectivenumberofparameters(pD),devianceinformationcriterion(DIC),
coefficientofdetermination(R2),adjustedcoefficientofdetermination(R2adj),relativestandarderror(Syx%)
andcorrectionfactor(CF)formodelsfitfromordinaryleastsquarewithlog-linearregressions.
tries, Inc, Boulder CO, USA). All codes used in this study parsimony, as indicated by the lower DIC values (Table 3).
werewrittenbytheauthors. TheinclusionoftheSGassignmentresultedinmodelswith
slightlylowerR2adjandhigherS comparedtothesame
yx%
modelstructurewithoutSG.
3 Results We observed differences with respect to the parameters
b and c among pioneer, mid- and late-successional species
3.1 Individualtreebiomassestimationmodelfits inmostofthemodelsthatincludedtheSGassignment(Ta-
bleS3andFig.S1).Thelate-successionalspeciestendedto
Although the NLS approach produced models with overall have higher intercepts and steeper slopes. Pioneer and mid-
higher values of R2 and R2adj and lower values of Syx%, successional species had lower differences in intercepts but
the DIC values indicated that the MOV and the OLS ap- stillstrongdifferencesintheslopes.
proachesproducedthebestmodels.ThemodelsM33(DBH Evaluations ofAGB predictionsfor individualtrees from
and WD as predictors) and M43 (DBH, SG and WD) were our two best models (as described in the Sect. 3.2) as well
the two best fitting models across all tree individuals (high fromthepantropicalmodel(Chaveetal.,2014)arepresented
R2 and R2adj and low Syx% and DIC values compared to in the Supplement of this study (Fig. S1). The models M33
other models). These two models also produced more reli- andM43hadlowerbiases(overestimationof0.6and3.5%,
able landscape predictions (see Sect. 3.2). The statistics for respectively)thanthetestedpantropicalmodel(underestima-
the goodness of fit for the 24 models are given in Table 3. tionof30%).
ForthemodelsfitwithOLS,whichrelyonlog-transformed
variables,theadditionofotherpredictorstogetherwithDBH 3.2 Landscape-levelbiomasspredictionsacross
systematicallydecreasedtheCFvalues.Thispatternsuggests scenarios
areductioninthebiasesresultingfrombacktransformation.
As expected, the addition of other predictors to a model To search for the model that best predicts AGB at the land-
containing only DBH systematically increased the models’ scapelevel,wetestedourmodels(excludingthosewithH as
Biogeosciences,13,1553–1570,2016 www.biogeosciences.net/13/1553/2016/
D.MagnaboscoMarraetal.:PredictingbiomassofcentralAmazonianforests 1561
Figure3.Predictedvs.observedabovegroundbiomass(AGB)alongsixforestscenarioscomposedof1001haplots.Thelineofequality
(1:1line)isshownasared/straightline.Forestscenariosweredesignedtoreflectlandscape-levelvariationsinfloristiccompositionand
size distribution of trees, typical of central Amazonian terra firme forests. Floristic composition and size-distribution scenarios followed
the sampling scheme described in Sect. 2.4.2 (Fig. 2) of this study. Models’ predictors: diameter at breast height (DBH) (cm), species’
successionalgroup(SG)(pioneers,mid-andlatesuccessional)andwooddensity(WD)(gcm−3).SeeTable2forthevariancemodeling
approachofdifferentequations.Notethatmodelscontainingtreetotalheight(H)aspredictorwereexcludedhere.
apredictor;Table2)acrossthe1001haplotsassembledfor We observed systematic biases ranging from −14%
eachofoursixforestscenarios(Figs.3–5)aswellasjointly (underestimation) to 38.8% (overestimation) in estimated
forallofthem(Fig.6). landscape-level AGB (Fig. 4). The models fit with NLS
The“true”AGBinour1haplots(fromthesummedmass tended to overestimate landscape-level AGB, with biases
of trees used to assemble the forest scenarios) varied from ranging from −3.6 up to 38.8%, both extreme values from
198.1 to 314.3 (early- to late-successional scenarios) and model series 1 (only DBH as predictor). Overall, the mod-
101.4 to 391.8Mgha−1 (small- to large-sized scenarios). els fit with NLS tended to capture changes in floristic com-
Theabilityofthevariousbiomassestimationmodelstopre- position better than in tree size distribution. The tested
dictthe“true”virtualbiomassvaluesgenerallyreflectedthe pantropicalmodelsystematicallyunderestimatedlandscape-
goodness of fit of the models for predicting individual tree level biomass, with a mean bias of −29.7% (Table S4 and
data(Table3andFigs.3–6).Thesamepatternwasobserved Fig.S2).
whenevaluatingthetestedpantropicalmodel,whichunder- The models fit with the OLS and particularly with the
estimatedboththeAGBofindividualtrees(Fig.S1)andin MOV approaches were clearly more efficient at capturing
allofourscenarios(TableS4andFig.S2). thevariationinfloristiccompositionandsizedistributionof
While some models produced accurate and satisfactory trees.Consequently,thesemodelsproducedthemostreliable
predictionsacrossallscenarios,otherssystematicallyunder- landscape-levelpredictionswithinthescenarios(Fig.3).As
oroverestimatedtheobservedAGB(Fig.3andFig.S2).The alsoindicatedbytheindividualtreemodelfits,theMOVap-
agreementbetweenmodelsandobservationswasinfluenced proach produced more reliable AGB predictions, especially
notonlybythedifferentcombinationsofpredictorsbutalso withmodelseries2and4.
bythedifferentmethodstomodelthevariance.Interestingly, In general, the models fit with the OLS and MOV ap-
despite producing the best fits to the individual tree data, proaches did not show systematic trends in under- or over-
models fit with NLS produced the least reliable landscape- estimation.ThemodelsfitwiththeOLSapproachhadbiases
level predictions, with model M11 (only DBH as predictor) ranging from −13.8 to 11.1%, with extreme values from
beingtheuniqueexceptionforthemid-andlate-successional model series 1 and 2, respectively. The models fit with the
scenarios(Fig.3). MOVapproachhadbiasesrangingfrom−14to10.5%,also
www.biogeosciences.net/13/1553/2016/ Biogeosciences,13,1553–1570,2016
1562 D.MagnaboscoMarraetal.:PredictingbiomassofcentralAmazonianforests
40 125
100
20 DBH 5705 DBH
0 ots 25
pl 0
ots40 ha 125
an bias (%) in 100 1−ha pl22400000 DBH + SGDBH + WD VariaNOMncLLOeSSV modeling −1RMSE (Mg ha) in 100 1−1112572570025055050050 DBH + SGDBH + WD VariaNOMncLLOeSSV modeling
Me n 0
40 a125
20 DBH + W Me17050 DBH + W
0 D + SG 2550 D + SG
0
Early−succ.Mid−succ.Late−succ.Small−sizedMid−sizedLarge−sized Early−succ.Mid−succ.Late−succS.mall−sizedMid−sizedLarge−sized
Floristic composition Size−distribution Floristic composition and size distribution
Figure5.Profilesrelatingtheroot-mean-squareerrorof12above-
Figure4.Profilesrelatingthebiasof12abovegroundtreebiomass
groundtreebiomassestimationmodelstestedalongsixforestsce-
estimation models tested along six forest scenarios composed of
narioscomposedof1001haplots.Forestscenariosweredesigned
1001haplots.Forestscenariosweredesignedtoreflectlandscape-
to reflect landscape-level variations in floristic composition and
level variations in floristic composition and size distribution of
sizedistributionoftrees,typicalofcentralAmazonianterrafirme
trees, typical of central Amazonian terra firme forests. Models’
forests.Models’predictors:diameteratbreastheight(DBH)(cm),
predictors:diameteratbreastheight(DBH)(cm),species’succes-
species’ successional group (SG) (pioneers, mid- and late succes-
sionalgroup(SG)(pioneers,mid-andlatesuccessional)andwood
density(WD)(gcm−3).Variancemodelingapproaches:nonlinear sional) and wood density (WD) (gcm−3). Variance modeling ap-
proaches:nonlinearleastsquare(NLS),ordinaryleastsquarewith
leastsquare(NLS),ordinaryleastsquarewithlog-linearregression
log-linear regression (OLS) and nonlinear with modeled variance
(OLS)andnonlinearwithmodeledvariance(MOV).Notethatmod-
(MOV).Notethatmodelscontainingtreetotalheight(H)aspre-
elscontainingtreetotalheight(H)aspredictorwereexcludedhere.
dictorwereexcludedhere.
withextremevaluesfrommodelseries1and2,respectively Asweexpected,theadditionofSGandWDimprovedthe
(Fig.4). qualityofthejointprediction.Thiswasevidencedbythesys-
The reported systematic biases led to strong differences tematicreductionofmodels’biasandRMSE.Notablyforthe
between the predicted and the observed AGB (Fig. 5). The NLSapproach,theinclusionofSGledtostrongreductionof
modelsfitwithNLSresultedinRMSEvaluesrangingfrom thebiasandRMSE(Fig.6).Interestingly,forthisapproach
16.8 up to 125.8Mgha−1. For the models fit with OLS, the addition of WD alone did not improve the estimations
the RMSE values ranged from 5.1 to 57.6Mgha−1. The accuracy.
MOVmodelshadRMSErangingfrom5.5to58.7Mgha−1.
The pantropical model’s predictions had a mean RMSE of
102.6Mgha−1(TableS4). 4 Discussion
By combining the bias and RMSE values, we could ob-
serve the overall models’ performance in predicting AGB 4.1 Individualtreebiomassestimationmodelfits
across scenarios (Fig. 6). When challenged to predict
biomass across all scenarios, the models fit with the MOV Thebest-performingallometrymodelstructuresforpredict-
approachproducedmorereliablepredictions(smallerrange ingthebiomassofindividualtreesincludedspecies-specific
ofbiasesandRMSE),exceptformodelseries1(onlyDBH predictors and either the OLS or MOV fitting approaches
asapredictor),forwhichtheOLSapproachperformedbet- (Fig. 3, Fig. 6 and Table S3). As we hypothesized, includ-
ter. Independently of applied predictors, the NLS approach ing both the SG and WD as predictors greatly increased
hadthehighestmeanandrangeofvaluesforbiasandRMSE. the models’ performance. When taken alone, adding either
Biogeosciences,13,1553–1570,2016 www.biogeosciences.net/13/1553/2016/