Developing a dynamic growth model for maritime pine in Asturias

25 Pages·2017·3.39 MB·English

Developing a dynamic growth model for maritime pine
in Asturias (NW Spain): comparison with nearby
regions
Manuel Arias-Rodil, Marcos Barrio-Anta, Ulises Diéguez-Aranda
To cite this version:
Manuel Arias-Rodil, Marcos Barrio-Anta, Ulises Diéguez-Aranda. Developing a dynamic growth
model for maritime pine in Asturias (NW Spain): comparison with nearby regions. Annals of Forest
Science, 2016, 73 (2), pp.297-320. ￿10.1007/s13595-015-0501-x￿. ￿hal-01356111￿
HAL Id: hal-01356111
https://hal.archives-ouvertes.fr/hal-01356111
Submitted on 24 Aug 2016
HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est
archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
AnnalsofForestScience(2016)73:297–320
DOI10.1007/s13595-015-0501-x
ORIGINAL PAPER
Developing a dynamic growth model for maritime
pine in Asturias (NW Spain): comparison
with nearby regions
ManuelArias-Rodil1 ·MarcosBarrio-Anta2·UlisesDie´guez-Aranda1
Received:29April2015/Accepted:28June2015/Publishedonline:25August2015
©INRAandSpringer-VerlagFrance2015
Abstract model may suffice for the whole of the NW Iberian
• Key message A dynamic growth model was developed Peninsula.
for maritime pine in Asturias. During the evaluation • Context Maritime pine is one of the most important
process, a stand volume ratio function proved the best tree species in NW Spain. There was no existing dynamic
oftwoalternativemethodsforestimatingmerchantable growthmodelforthisspeciesinAsturias.
volume.Comparisonofthedevelopedmodelwithexist- • Aims To develop a dynamic growth model for maritime
ing models for nearby regions showed that a single pineinAsturias,byevaluatingtwodifferentmethodsofesti-
matingvolume(adisaggregationsystemandastandvolume
ratio function), and to compare the developed model with
existing models for Galicia and northern Portugal are the
Communicatedby:AaronRWEISKITTEL
goalsofthisstudy.
HandlingEditor:AaronRWEISKITTEL • Methods The dynamic model is based on the state-space
approach, in which three state variables characterize the
Contribution of the co-authors Manuel Arias-Rodil: did the
analyses,implementedthestandgrowthsimulatorinR,andwrote stand at any point in time: dominant height, number of
themanuscript. stemsperhectareandstandbasalarea.Thetransitionfunc-
MarcosBarrio-Anta:advisedindiscussionandrevisedthetext.
tion for the first variable was developed on the basis of
Ulises Die´guez-Aranda: supervised the work, advised in discus-
stem analysis data in a previous study, while the corre-
sion,helpedintheimplementationofthestandgrowthsimulator
inR,andrevisedthetext sponding functions for the last two variables were simu-
ltaneously fitted with data obtained from successive mea-
(cid:2) ManuelArias-Rodil surements of permanent plots. An appendix outlining the
[email protected] implementation of a stand growth simulator in the R
MarcosBarrio-Anta environmentisincludedtofacilitatemodeluseandevalua-
[email protected] tion.
• Results When the whole model was used to project the
UlisesDie´guez-Aranda
standconditions,thestandvolumeratiofunctionperformed
[email protected]
best, yielding a root mean square error of 22.4 m3 ha−1
and a critical error of 18.4 %. Comparison with models
1 UnidaddeGestio´nForestalSostenible(UXFS), developed for other regions revealed both similarities and
DepartamentodeIngenier´ıaAgroforestal,Universidad
differences,someofwhichmaybeattributedtoanunequal
de Santiago de Compostela, Escuela Polite´cnica Superior,
C/BenignoLedo,CampusUniversitario,27002Lugo,Spain distribution of the available data in age and site quality
classes.
2 Grupo de Investigacio´n en Sistemas Forestales Atla´nticos •ConclusionTheproposeddynamicgrowthmodelprovided
(GIS-Forest),DepartamentodeBiolog´ıadeOrganismos
accurateresults,andcomparisonwithotherregion-specific
ySistemas,EscuelaPolite´cnicadeMieres,Universidad
deOviedo,C/GonzaloGutie´rrezQuiro´ss/n, modelsshowedthatasingledynamicmodelmaysufficefor
Mieres33600,Spain thewholeoftheNWIberianPeninsula.
298 M.Arias-Rodiletal.
Keywords PinuspinasterAit.·State-spaceapproach· single-species plantations, whole-stand models therefore
Disaggregationsystem·Growthsimulator·Whole-stand represent a good compromise between accuracy and gen-
model erality (Garc´ıa 2003). They use stand variables as inputs
and outputs, although some of the models may include
a mathematical disaggregation of the stand-level inputs
1Introduction to estimate the number of stems in different diameter
classes (referred to as diameter distribution models by
The growing stock of maritime pine (Pinus pinaster Ait.) BurkhartandTome´ (2012),p.234).
is the largest in Spain, currently representing 15 % of Volumepredictiontoanymerchantablelimit(fordiffer-
the timber volume and 27 % of the annual harvested ent timber assortments) is commonly achieved (Burkhart
volume (MAGRAMA 2010). The wood is commonly and Tome´ 2012, p. 270) by using (i) a stand volume ratio
used to produce sawn timber and also in the pulp or function,whichprovidesthevolumetoatopdiameterlimit
wood-based panel industries, depending on log dimen- (Barrio-Anta et al. 2008) or to both a top diameter limit
sions (Sanz et al. 2006). Maritime pine covers a total and a diameter class (Amateis et al. 1986), or (ii) a dis-
area of 22,500 ha in the region of Asturias (NW Spain) aggregation system, based on the joint use of a diameter
(5%,MAGRAMA2012). distribution, a height–diameter (h–d) relationship, and a
Because of the importance of this species in Asturias, stem taper function, to estimate the merchantable volume
several studies have addressed the following aspects: site accordingtotopdiameterlimitsandloglengths.Thelatter
quality (A´lvarez-A´lvarez et al. 2011; Arias-Rodil et al. approach has been used in many models (Die´guez-Aranda
2015a), edaphic factors and nutritional status (Afif Khouri et al. 2006; Castedo-Dorado et al. 2007; Go´mez-Garc´ıa
et al. 2009), tree biomass (Canga L´ıbano et al. 2009) et al. 2014a), and it is the approach most commonly used
and diameter distribution (Gorgoso Varela et al. 2009). to estimate volume within a dynamic growth model. Nev-
However, no growth and yield models have yet been ertheless, there are several disadvantages associated with
developed for the species. Such models are very use- thismethod,suchaspossiblecomputationalinefficiency(it
ful for forest practitioners as they enable stand growth depends on several functions, some of which may involve
simulation for different management practices. This type iterative procedures) and the inclusion of several submod-
of information is essential for decision-making pur- els, with their corresponding estimation errors. Therefore,
poses in the context of sustainable forest management it is useful to evaluate whether a stand volume equation
(Vanclay 1994, p.xv). Empirical models are widely used can yield results that are as accurate as those produced by
as practical tools in forest management, as they require the disaggregation system, but in a more computationally
few and easy-to-measure (inexpensive) input variables and efficientway.
provide accurate yield predictions. However, as they are Growth and yield models have already been devel-
not linked to the underlying growth processes, they do not oped for maritime pine plantations in nearby Atlantic
takeintoaccountsuddendisturbances(e.g.extremeweather regions: the models developed for the species in Gali-
events) or long-term environmental changes (i.e. climate cia are summarized in Die´guez-Aranda et al. (2009)
change). A detailed review of model types can be found and for northern Portugal (Taˆmega Valley, model termed
inMa¨kela¨etal.(2000),Monserud(2003)andPretzschetal. ModisPinaster) in Fonseca (2004) and Fonseca et al.
(2008). (2012). These models provide information about the
Empirical growth models can be classified into three behaviour of this species in each region, enabling com-
types (Davis et al. 2001, p. 185): (i) whole-stand mod- parison with behaviour of the species in Asturias. We
els, (ii) size-class models and (iii) individual-tree mod- consider such comparison of interest, as it represents
els. The differences between these models are based another way of validating the proposed dynamic growth
on the level of resolution of predictions and data model.
requirements. Model types (ii) and (iii) are useful as Themainobjectivesofthisstudywereasfollows:(i)to
research tools (they describe stand behaviour), but over- develop a dynamic growth model for maritime pine in
parameterization is common and limits the accuracy and Asturias and to evaluate its performance in total and mer-
precisionofpredictionsatforestmanagementlevel(Garc´ıa chantable stand volume estimation, (ii) to compare the
2003). In addition, collection of individual-tree and size- predictionoftotalandmerchantablestandvolumebyusing
class data is more expensive than collection of data at both a stand volume ratio equation and a disaggrega-
stand level. In the special case of the management of tion system and (iii) to compare the performance of the
DynamicgrowthmodelformaritimepineinAsturias 299
proposed dynamic growth model with those of existing description; observed and predicted trajectories in Fig. 1,
models for Galicia and northern Portugal, as regards pro- topleft),whichwasthereforeusedasthetransitionfunction
jection of stand variables, estimation of total stand volume fordominantheightandforcomputingthesiteindexofeach
and determination of optimal biological rotation age (the plot-inventorycombination.ItsexpressionisshowninEq.1
age of maximum mean annual increment—MAI—of total andthecorrespondingparametersareshowninTable2.Site
standvolume). indexvalueswerethenaveragedbyplot,assumingthatsite
qualityisconstantovertime.
2Materialandmethods H = (cid:2) b1 (cid:3)(cid:2) (cid:3) (1)
2 b t b3
1− 1− 1 1
2.1Data H t
1 2
Thedatausedtodevelopthemodelwereobtainedfromtwo whereH isdominantheight(m)ataget (years),obtained
2 2
networksofplotsinstalledinpure,even-agedmaritimepine from dominant height H at age t , and b and b are
1 1 1 3
stands: (i) 74 permanent plots and (ii) 18 plots of a thin- theparametersoftheADAformulationofHossfeld(1822)
ningtrial.Thepermanentplotsofthefirstdatasourcewere model,presentedbyMcDillandAmateis(1992).
establishedandmeasuredin2007.Theplotswereinstalled Transition functions serve to describe the natural evolu-
throughouttheareaofdistributionofthespeciesinAsturias tionofstands,andthereforegrowthintervalsusedforfitting
(mainly in the NW of the region), covering the existing should not include silvicultural treatments (e.g. thinning)
range of ages, stand densities and site qualities. The plot or random disturbances (e.g. fire, wind damage). In accor-
size ranged from 700 to 900 m2 to include a minimum of dance with this criterion, 37 growth intervals in the first
30trees.Asecondmeasurementwasmadeonasubsetof58 data source and 8 in the second (4 between the first and
of the 74 initially established plots (in 2011 and 2012), as the second inventories and 4 between the second and the
somedisappearedasaresultofforestfiresorclear-cutting. third)weredisregarded.Thus,thetransitionfunctionswere
The18plotscorrespondingtotheseconddatasourceare finallyfittedfor88plot-inventorycombinations(35,39and
locatedin6sites(3plotsof1000m2persite),inwhicheach 14 plots for the first, second and third inventories, respec-
plotwassubjectedtoadifferentthinningtreatment(nothinning, tively).Inaddition,allthe186plot-inventorycombinations
control; light low thinning and heavy low thinning). Three available (92, 76 and 18 for the first, second and third
measurements were carried out (in 2009, 2011 and 2013), inventories, respectively) were used to develop part of the
whichimpliestwoavailablegrowthintervalsforeachplot. disaggregation system (diameter distribution and height–
In each plot corresponding to the first and second data diameterrelationship)andthestandvolumeratiofunction.
sources, diameter at breast height (d, cm, at 1.3 m from Summary statistics of the stand variables used in model
the ground) and total height (h, m) were measured to the developmentaregiveninTable1.
nearest0.1cmand0.1m,respectively,inalltrees.Descrip-
tive variables were also recorded for each tree, e.g. if they 2.2Dynamicmodel
werealiveordead.Thefollowingstandvariableswerecal-
culated: age (t, years), dominant height (H, m, defined as Thedynamicwhole-standmodeldevelopedinthisstudyis
the mean height of the dominant trees, i.e. the 100 thick- based on the state-space approach (Garc´ıa 1994), which is
esttreesperhectare),siteindex(S,definedasthedominant suitable for systems that evolve over time, such as forest
height of the stand at a reference age of 20 years), num- stands.Thecurrentstateofastandcanbedefinedbyalist
berofstemsperhectare(N),standbasalarea(G,m2ha−1), of state variables, which can be projected to the future by
and total (V, m3 ha−1) and merchantable (V ) volumes to thetransitionfunctions.Theideaistocharacterizethestate
i
differenttopdiameterlimits(from1to35cmby1-cminter- of the system at any point in time so that the future state
vals). The stand volumes were obtained by aggregation of doesnotdependonthepaststate.Threestatevariablesare
the corresponding tree volumes estimated using the stem often used (dominant height, number of stems per hectare
taperfunctionfittedbyArias-Rodiletal.(2015b).Inaddi- and stand basal area, e.g. Die´guez-Aranda et al. 2006;
tion, a dominant height projection function was developed Castedo-Doradoetal.2007;A´lvarez-Gonza´lezetal.2010;
byA´lvarez-A´lvarezetal.(2011)onthebasisofstemanaly- Go´mez-Garc´ıa et al. 2014a), although for high-intensity
sisdatafrom146treeslocatedclosetotheplotsofthefirst treatments (e.g. high-height pruning, heavy thinning), use
data source (see A´lvarez-A´lvarez et al. (2011) for the data of a fourth variable, such as a measure of stand closure, is
300 M.Arias-Rodiletal.
Fig.1 Observedtrajectories Observation Prediction Asturias
(solidgreylines)ofdominant Region−specific prediction Prediction Portugal
height(H,left)andstandbasal 20 40 60
area(G,right)byregion,and H G
thecorrespondingheight 30 100
predictions(siteindicesof7,11,
15and19m)andstandbasal 25 80
area(valuesof10,25,40and
55m2ha−1at20years), 20
providedbytheregion-specific s 60s
a a
model(solidblacklines),the 15 uri uri
submodelsforAsturias(dashed Ast 40Ast
blacklines),andthesubmodels 10
forPortugal(dottedblacklines).
20
ThePortuguesestandbasalarea 5
submodeldependsonstem
densityatbothinitialand 0 0
projectionage,aswellason 30 H 100 G
standbasalareaattheinitial
age;thegraphwasobtained 25
80
fpartaos2ms0estyheteharrmosuo(gtrhthael1imt0y0oc0rtuasrltvietemytcshuahrtav−e1 −21m ha) 20 st) 60st)
ahtthhlesreiooghuhdegtei,hgpwhe1thn1igdcmrshoowawnttah2ds0ocomuybreitvnaaeairnsnp)etadssfirnogm ue (H, m; G, 1105 Galicia (coa 40Galicia(coa
Val 5 20
0 0
H G
30 100
25
80
20 erior) 60erior)
15 a (int a(int
ci 40ci
10 ali ali
G G
20
5
0 0
20 40 60
Age (years)
recommended to account for site occupancy and to take of projecting from t to t and then from t to t must be
1 2 2 3
into account the response to these treatments (e.g. Garc´ıa the same as projecting directly from t to t ) and causal-
1 3
2011, 2013, Garc´ıa et al. 2011). In this study, we consid- ity (a change in the state can only be influenced by inputs
eredsingle-speciesstands,moderatelythinnedfrombelow, within the relevant time interval). Once the state variables
andundertheseconditions,H,N andGareabletoprovide areknownforagiventime,totalandmerchantablevolumes
agooddescriptionofthestandconditionatanyage.Three canbeestimatedindifferentways.Onewayistouseastand
transition functions enable projection of the state variables volume ratio function (e.g. Tewari et al. 2014) to estimate
overtime.AccordingtoGarc´ıa(1994),transitionfunctions totalandmerchantablevolumedirectlyfromthestatevari-
must possess certain desirable properties: consistency (no ables.Analternativeapproachistouseadisaggregationsys-
change for zero elapsed time), path-invariance (the result tem(e.g.Die´guez-Arandaetal.2006,Castedo-Doradoetal.
DynamicgrowthmodelformaritimepineinAsturias 301
Table1 Datasummary
Variable Mean Min. Max. Std.dev. Mean Min. Max. Std.dev. Mean Min. Max. Std.dev.
Firstdatasource Firstinventory(74plots) Secondinventory(58plots)
t 29.9 8.0 61.5 14.7 32.1 12.0 63.0 13.0
H 15.5 5.5 27.0 5.7 16.9 9.8 29.3 4.8
N 1065 378 2480 469 943 111 1900 421
G 38.6 7.8 76.2 16.9 39.8 12.5 75.0 14.1
V 263.0 23.5 742.4 179.9 278.2 80.9 785.4 156.5
S 11.9 7.1 19.4 2.2 11.9 7.1 19.4 2.3
Seconddatasource Firstinventory(18plots) Secondinventory(18plots) Thirdinventory(18plots)
t 17.5 12.0 29.0 6.5 19.5 14.0 31.0 6.5 21.5 16.0 33.0 6.5
H 10.6 7.8 15.7 3.0 12.0 9.2 17.0 2.9 13.3 10.0 18.3 2.9
N 993 470 1490 341 984 470 1450 336 908 460 1430 312
G 20.6 13.7 40.6 7.1 24.6 16.6 44.1 7.2 27.6 21.0 42.0 6.0
V 100.0 53.7 245.6 52.8 127.1 75.3 278.0 53.0 154.2 97.0 278.2 48.6
S 12.7 11.0 15.1 1.3 12.7 11.0 15.1 1.3 12.7 11.0 15.1 1.3
tstandage(years),H dominantheight(m),Nnumberofstemsperhectare,Gstandbasalarea(m2ha−1),V totalstandvolume(m3ha−1),Ssite
index(m,atreferenceageof20years)
2007, Go´mez-Garc´ıa et al. 2014a): diameter distribu- bestmodelfordescribingmortalitywasthatofTome´ etal.
tion models generate the number of stems in diameter (1997):
classes, while generalized height–diameter models predict
theheightfortheaveragetreeofeachclass,andtaperfunc-
N =N exp(a (t −t )) (2)
2 1 0 2 1
tionsareusedtocomputeandclassifythevolumeaccording
to top diameter limits and log lengths, which are specified
bymarketrequirements. where N2 is the projected number of stems per hectare
In the following sections, we explain how we devel- at age t2 (years), N1 and t1 are respectively the initial
oped the transition functions, the stand volume ratio values of number of stems and age and a0 is the model
function and the disaggregation system. As already men- parameter.
tioned, A´lvarez-A´lvarez et al. (2011) have already devel- Use of the ADA and GADA approaches implies that
oped the height growth curves using stem analysis data for prediction, the initial value of the variable must be
(1), and therefore, these were not re-fitted in this study. known at a given age. This value is commonly known for
Most of the submodels were fitted by the ordinary least the number of stems (which are easily counted), but not
squares (OLS) technique, by using the nls function of R for the stand basal area. If this variable is not known, an
(RCoreTeam2015). initialization equation is used to predict stand basal area
fromotherstandvariables.Therefore,wedevelopedastand
2.2.1Mortalityandstandbasalareagrowthfunctions basal area growth system comprising two functions: one
for projection and another for initialization. For modelling
The algebraic difference approach (ADA—Bailey and standbasalareagrowth,previousanalysesshowedthatthe
Clutter 1974) and its generalization (GADA—Cieszewski bestmodelwastheGADAformulationofHossfeld(1822)
andBailey2000)wereusedtodevelopthetransitionfunc- model (3). For the initialization function, previous analy-
tion for number of stems per hectare and stand basal area. ses showed that basal area is allometrically related to S,
Thetransitionfunctionforthedecreaseinnumberofstems N and t (4). In addition, the model was fitted with data
considers only the natural mortality (i.e. that caused by from plots younger than 20 years, given that its predic-
competition for light, water and soil nutrients—Peet and tive capacity was found to be much lower when the whole
Christensen 1987). Preliminary analysis showed that the data set was used. Thus, for stands older than 20 years,
302 M.Arias-Rodiletal.
we assumed that stand basal area should be obtained from quadratic mean diameter and a top diameter limit. Both
inventorydata. werecombinedinthesameequation,giventhattotalstand
volume becomes a special case of the stand volume ratio
X
G = 0 (3) equation when the top diameter limit is equal to zero
2 1+ b2t−b3 (GregoireandSchabenberger1996).Thefinalequationhas
X 2
0 thefollowingform:
(cid:2) (cid:4) (cid:3) (cid:5) (cid:6)
where X0 = 21 G1+ G21+4b2G1t1−b3 ; G2 is the Vi =VRi =a0Ga1Ha2exp b0Dgb1dib2 (5)
projected stand basal area (m2 ha−1) at aget 2 (years); where Vi is the stand volume (m3 ha−1) to a top diame-
G1 and t1 the initial values of stand basal area and age, ter limit di (cm), G the stand basal area (m2 ha−1), H the
respectively; and b2 and b3 the model parameters (b1 was dominantheight(m),Dg thequadraticmeandiameter(cm)
non-significant). and a -a and b -b are the model parameters. Note that
0 2 0 2
V = a0Ga1H(cid:5)a2 represen(cid:6)tsthetotalstandvolumeequation
G1 =a0Sa1Na2ta3 (4) andR =exp b Db1db2 theratioequation.
i 0 g i
Intheprocessoffittingbothtransitionfunctions,weused
the base-age invariant dummy variables method proposed Disaggregation system After selecting an appropriate
by Cieszewski et al. (2000), which estimates site-specific functionalformtocharacterizethediameterdistribution,we
effectsundertheassumptionthatdatameasurementsalways obtain the parameter estimates and predict the number of
include measurement and environmental errors (both on stems per diameter class. The height of the average tree in
the left- and right-hand sides of the model) that must be eachdiameterclassispredictedusingtheh–d relationship.
modelled. Previous fitting showed that the residuals of Total and merchantable tree volume as dependent on top
both the transition function for stand basal area and num- diameter and log length are predicted using a stem taper
ber of stems were correlated (Pearson’s moment product function and aggregating the individual-tree results to the
correlation of 0.24). To account for this correlation, we wholestand.
fitted the models simultaneously by using the seemingly
unrelated regression (SUR—Zellner 1962) technique for
Diameter distribution We used a two-parameter Weibull
nonlinear models. This method is completed in two steps,
function to model the diameter distribution, given that it is
which are based on the OLS technique: (i) separate fit-
the most frequently used in forest growth models because
ting of each of the equations considered in the system
ofitsflexibilityandsimplicity(Maltamo1995;Kangasand
(the two transition functions in this case) and (ii) re-fitting
Maltamo2000;Torres-Rojoetal.2000).Theexpressionof
simultaneously all equations of the system considering
thetwo-parameterWeibulldensityfunctionisasfollows:
the correlation between residuals obtained in the first (cid:5) (cid:6)(cid:5) (cid:6) (cid:5) (cid:5) (cid:6) (cid:6)
c x c−1 x c
step. f(x)= exp − (6)
The first step of SUR fitting was accomplished by the b b b
OLStechnique,whileforthesecondstep,theresidualsum- wherex istherandomvariable,andbandc(bothpositive)
of-squares computation of the system was implemented in arethescaleandshape(skewness)parameters,respectively.
afunctionandthenminimizedbythenlmfunctionofR(R The parameters of the Weibull function can be obtained
CoreTeam2015). bydifferentmethods,whichcanbeclassifiedintwogroups:
parameterpredictionandparameterrecovery(Hyink1980;
2.2.2Volumeestimation Vanclay 1994 p. 23). Considering the scope of this study,
parameter recovery was used as it proved the best method
Standvolumeratiofunction Predictingthemerchantable in several previous studies (Cao et al. 1982; Reynolds
volume (to a top diameter limit) as a ratio of total vol- et al. 1988; Torres-Rojo et al. 2000). The parameters are
ume was originally proposed by Burkhart (1977) at tree directlyobtainedfrompercentiles(CaoandBurkhart1984)
level, for which a total volume and a ratio equation must or moments (Newby 1980; Burk and Newberry 1984) of
be fitted. Barrio-Anta et al. (2008) applied this approach the diameter distribution, estimated from its relation with
to the estimation of merchantable stand volume. As in stand variables. Recovering parameters from the moments
the latter study, we considered a total stand volume equa- (known as the moments method) is the only method that
tionbasedonanallometricrelationshiptostandbasalarea directlyensuresthatthesumofthedisaggregatedbasalarea
and dominant height, and a ratio equation depending on obtainedbytheWeibullfunctionequalsthestandbasalarea
DynamicgrowthmodelformaritimepineinAsturias 303
provided by an explicit growth function of this variable, (Eq. 10), which was also used by Go´mez-Garc´ıa et al.
resultinginnumericcompatibility(Hyink1980;Kangasand (2014b).
(cid:2) (cid:2) (cid:3)(cid:3)
Maltamo2000;Torres-Rojoetal.2000);wethereforechose
1 1
tousethismethod. h=1.3+(H−1.3)exp (a0+a1H +a2Dg) −
d D
d
Thefunctionparameterscanberecoveredfromthefirst
(10)
two moments of the diameter distribution (mean and vari-
ance)(Caoetal.1982).Thearithmeticmeandiameter(D ,
m wherehisthetreeheight(m),d thetreediameteratbreast
cm) corresponds to mean, and the variance (var) is esti-
height (cm), H the dominant height (m), D the dominant
mated as var = D2 − D2 (D , quadratic mean diameter, d
g m g diameter (cm), D the quadratic mean diameter (cm) and
g
cm).TheknownvaluesofD andD canbeusedtoobtain
g m a -a aremodelparameters.
0 2
var and Eq. 7 can be solved iteratively for parameter c.
Practical use of the generalized h–d equation requires
Finally,parameterb iscomputeddirectlyfromEq.8using
estimation of the dominant diameter, which was estimated
previouslyobtainedvaluesofD andtheparameterc.
m fromthediameterdistribution,asitisdifficulttoprojectin
(cid:7) (cid:2) (cid:3) (cid:2) (cid:3)(cid:8)
D2 2 1 time(Lappi1997).
var= m (cid:2) 1+ −(cid:2)2 1+ (7)
(cid:2)2(1+1/c) c c
Stem taper function When the height of the average tree
D in each diameter class is known, the merchantable volume
b= m (8)
(cid:2)(1+1/c) fromthestumptoafixedtopdiameterlimitcanbeestimated
directlybyusingatreevolumeratioequationorbyintegra-
where D the arithmetic mean diameter (cm), b and c the
m tion of a stem taper function (Castedo Dorado et al. 2005;
Weibullparametersand(cid:2)theGammafunction.
Die´guez-Aranda et al. 2006). According to these authors,
The quadratic mea(cid:9)n diameter(cid:10)(Dg)(cid:11)can be obtained thelatteristhemostcommonlyusedapproach.
directlyfromNandG D =100 4G ,butD remains In this study, we used the stem taper function fitted in
g m
πN Arias-Rodil et al. (2015b), which is based on the Kozak
unknown. However, it can be modelled through a rela- (2004) model. Assuming that no additional measurements
tionship with the quadratic mean diameter and other stand will be available to calibrate the model, we considered
variables. The best set of stand variables were included in the parameter estimates corresponding to the fixed-effects
the followingexpression, which ensures thatDm is always model,whichwasfittedbyOLS:
lowerthanorequaltoD (Frazier1981):
g
Dm =Dg−exp(a0+a1H +a2N) (9) di =a0da1ha2x(b1(hi/h)4+b2(1/exp(d/h))+b3x0.1+b4(1/d)+b5hw+b6x)
(11)
where D is the arithmetic mean diameter (cm), D the
m g
quadraticmeandiameter(cm),H thedominantheight(m),
N thenumberofstemsperhectareanda -a arethemodel whered (cm)isthetopdiameteratstemheighth (m),dthe
0 2 i i
parameters.ThisequationwasfittedbytheOLStechnique. diameter at breast height (cm), h the total tree height (m),
w
w = 1 − (h/h)1/3, x = , anda -a and
i 1−(1.3/h)1/3 0 2
Height–diameter relationship A height–diameter rela- b1-b6arethemodelparameters.
tionship was used to estimate the tree height for each Equation 11 cannot be analytically solved for h and
i
diameter class. A single h–d relationship may not be ade- cannot be directly integrated to obtain the volume to
quateforallsituationsbecauseitvariesbetweenstandsand a top diameter limit (v), which implies that numerical
i
alsowithage(Curtis1967).Tosolvethis,thegeneralizedh– procedures should be used. The h at which a specific
i
drelationshipsincludestandvariablestolocalizetheheight diameter (d) is reached was obtained by the optimize
i
predictionsforeachstandage. function of R (R Core Team 2015), which uses the
Inmodellingtheh–drelationship,weselectedequations method of Brent (1973), considered appropriate for one-
that are constrained to predict stand dominant height (H) dimensional optimization. On the other hand, the total and
fromdominantdiameter(D )andtoyield1.3mforzerod merchantable tree volumes (v and v, respectively) of the
d i
values (to prevent negative estimates in small trees). From average tree of each diameter class were estimated by
these, previous analyses revealed better performance of a numericalintegrationwiththeintegratefunctionofthe
generalizedformofthemodelofBurkhartandStrub(1974) same software. The products of these tree volumes and
304 M.Arias-Rodiletal.
the corresponding number of stems were finally aggre- is generally considered realistic and reasonable in forest
gated to estimate total (V) and merchantable (V) stand growthmodelling(Huangetal.2003).
i
volumes. (cid:12)
(cid:13)
(cid:13) (cid:15)n
2.3Overallevaluationofthemodel τ(cid:14) 1 (y −yˆ )2
χ2 i i
crit i=1
E = (14)
Eachsubmodelofthedynamicmodelwasevaluatedgraph- crit y¯
ically and numerically. We visually inspected plots of
whereτ isastandardnormaldeviateatthespecifiedprob-
residuals against estimated values, to evaluate the pres-
ability level (τ = 1.96 for α = 0.05), χ2 is the value of
ence of heteroscedasticity, and graphs of the predictions crit
the χ2 distribution obtained for α = 0.05 and n degrees
of transition functions Eqs. 1, 2 and 3 overlaid on the
of freedom, and the remaining variables are as previously
observed trajectories, to assess the adequacy of each sub-
defined.
model. Numerical analyses were based on the root mean
Within the disaggregation system, we applied the
square error (RMSE) and the adjusted coefficient of deter-
minationfornonlinearmodels(R2 ).Bothofthesetakeinto Kolmogorov–Smirnov test (KS): (i) to all plot-inventory
adj combinations,toassessthesuitabilityofthetwo-parameter
accountthenumberofparameters,thuspenalizinginflated
Weibull function for describing the diameter distribution
modelparameterization.Theseerrorstatisticsareexpressed
and (ii) to all projected stands (projected as explained in
asfollows:
thepreviousparagraph),toevaluatetheperformanceofthe
(cid:12)
(cid:13)(cid:13)(cid:15)n parameter recovery approach when projecting state vari-
(cid:13)(cid:13) (yi −yˆi)2 ables. The KS test compares the estimated and the real
RMSE=(cid:14)i=1 (12) diameter distribution. Because the estimated distribution
n−p parameters are determined from the data, Lilliefors (1967)
statedthattheKS-statisticexistingdistributionisnolonger
valid and should be obtained by Monte Carlo simulation.
(cid:15)n For each plot-inventory combination or projected stand,
(y −yˆ )2 we generated 10,000 independent identically distributed
i i
R2 =1− i=1 n−1 (13) pseudo-randomsamplesunderthenullhypothesis(i.e.with
adj (cid:15)n n−p recovered parameters), and we computed the KS statistic
(yi −y¯)2 for each sample. This enables approximation of the KS-
i=1 statistic distribution under the null hypothesis, which was
subsequently used: if the KS statistic value obtained from
where yi and yˆi represent the observed and estimated val- thecomparisonbetweentheestimatedandrealdistribution
ues, y¯ is the average of yi values, n is the total number of of a plot exceeds the critical value at a specified signifi-
observations,andpisthenumberofmodelparameters. cancelevel(obtainedfromtheapproximateddistributionof
Althougheachsubmodelwouldperformadequately,this KSstatistic),thehypothesisthattheobservationsbelongto
does not guarantee that the overall dynamic growth model aWeibulldistributionofthespecifiedparametersshouldbe
providesreliableresults.Thus,anoverallevaluationshould rejected.
bemade.Becauseofthelackofanindependentdataset,we
projected the information (H, N and G) of the first inven- 2.4SimulatorinR
tory(35plots)totheagesofthesecondandthirdinventories
byusingthetransitionfunctions.Attheseages,theprojec- The dynamic growth model of Asturias can be included
tionswereusedtoestimatethetotalandmerchantablestand in a stand growth simulator, which enables simulation of
volumes,astheyarecloselyrelatedtotheeconomicvalueof different management schedules (depending on the tim-
astand.Forthispurpose,weappliedboththestandvolume ing, intensity and type of cutting). Simulation of thinning
ratioapproachandthedisaggregationsystem. operationsvaries dependingon whether the disaggregation
Thetotalstandvolumeestimates andstatevariablepro- system (information about diameter classes) or the stand
jections were also evaluated in terms of the critical error volume ratio function is used. For implementing this sim-
(Reynolds 1984), which was computed by rearranging the ulator, the R statistical software (R Core Team 2015) is
statistic of Freese (1960) (14), and are expressed as a per- a good means of transfer for the following reasons: (i) it
centageoftheobservedmean.Acriticalerrorof10to20% is commonly used by statistical analysts, (ii) it enables a
DynamicgrowthmodelformaritimepineinAsturias 305
testing workspace and (iii) it is easy to learn to use. Rea- (500,900,1300and1700stemsha−1),toobtain16exam-
sons (i) and (ii) ease the adaptation and expansion of the ple stands. To use the same initial stands for comparison,
code by the research community, while reason (iii) facili- initial basal area was computed from these variables and
tates the use of the dynamic growth model by forest prac- from t equal to 20 years, using Eq. 4 (with parameters to
titioners. More details about thinning simulation and code be obtained in the fitting step). Region-specific dynamic
structureareshownintheAppendix. growth models were then used to simulate the stand evo-
lution following a no-thinning schedule. Finally, for each
2.5Comparisonwithotherdynamicmodels stand and model, the age of maximum total stand volume
MAIwasobtained,i.e.thetimeatwhichthebiologicalpro-
Once the dynamic growth model was developed for mar- ductivity is highest or the optimal biological rotation age
itimepineinAsturias,itwascomparedwithmodelsprevi- (AveryandBurkhart2002,pp.353–355).
ously developed for the same species in Galicia (included
in Die´guez-Aranda et al. 2009) and northern Portugal
(ModisPinaster,Fonseca2004;Fonsecaetal.2012).These 3Results
modelsarebasedonthesamestatevariablesasthepresent
model and include a disaggregation system (they do not 3.1Dynamicmodel
have a stand volume ratio equation). The Galician model
includes dummy variables in the dominant height and Table 2 shows the parameter estimates obtained during
stand basal area transition functions, as well as in the h–d the fitting step for all the submodels included in the pro-
relationship, to differentiate between two ecoregions (Al´ıa poseddynamicgrowthmodelformaritimepineinAsturias.
Miranda et al. 2009): coast and interior. No mortality was Note that both height growth and stem taper functions
observed in Pinus pinaster stands in Galicia; thus, the Eqs. 1 and 11, respectively, have already been developed
model does not include a stem density reduction equation. (no error statistics shown in Table 2). For transition func-
ModisPinaster uses a disaggregation system based on the tions, height and stand basal area growth functions are
Johnson’s S distribution, a h–d relationship, and a tree shown in Fig.1 (top left and top right, respectively),
B
volume ratio function to estimate total and merchantable while mortality curves are shown in Fig. 2 (solid grey and
volumes.Weusedanewlydevelopedh–d relationshippro- black lines). The predicted trajectories of these transition
posed by Go´mez-Garc´ıa et al. (2015) for this species and functions were plotted for different initial conditions (val-
region,insteadofthatconsideredbyFonseca(2004). ues of H, N and G) at age 20 years, overlaid on the
Themodelcomparisonbetween regionswas carried out observations.
onthebasisof(i)projectionofstatevariables(H,NandG),
(ii) prediction of diameter distribution and total stand vol- 3.2Overallevaluationofthemodel
umefromprojectedstatevariablesand(iii)theageatwhich
the mean annual increment (MAI) of total stand volume is Table3showstheRMSEandcriticalerrorsobtainedwhen
maximal(optimalbiologicalrotationage).Forcomparisons using the developed dynamic growth model for projection
(i)and(ii),weappliedeachregion-specificdynamicgrowth ofstatevariablesandestimationoftotalstandvolume.The
model to the first-inventory plots of the data set used in latterwascomputedusingboththestandvolumeratiofunc-
thepresentstudy,projectingtheirstatevariablestotheages tion and the disaggregation system. The results of apply-
of second and third inventories, and then predicting diam- ing other region-specific dynamic models to our data are
eter distribution and estimating the total stand volume, as also included but will be discussed later. Additionally, for
donefortheoverallevaluationofthedynamicgrowthmodel comparisonbetweenvolumeestimationalternatives,Fig.3
for Asturias. We subsequently calculated the RMSE and showstheRMSEvaluesinmerchantablevolumeestimation
critical errors in projecting the state variables and in esti- for0-to30-cmtopdiameterlimits,obtainedbyusingboth
matingtotalstandvolume.Forthediameterdistribution,we approaches. The stand volume ratio function performed
computedthemeanandvarianceofplot-levelKolmogorov– slightly better and was therefore used in comparison with
Smirnov statistics, obtained from the comparison between otherdynamicmodels.
predicted(byeachregion-specificmodel)andrealdiameter The Kolmogorov–Smirnov test showed that the Weibull
distribution. functionsuccessfullyexplained(ata5%significancelevel)
For comparison (iii), we combined four site indices (7, the diameter distribution of Pinus pinaster in 94.1 % of
11, 15 and 19 m) and four stem densities at 20 years the plot-inventory combinations, and when evaluating the