Base-age invariant models for predicting individual tree accumulated annual resin yield using two tapping methods in maritime pine (Pinus pinaster Ait.) forests in north-western Spain

Fo es Ecology and Managemen 549 (2023) 121501 A ailable online 19 Oc obe 2023 0378-1127/© 2023 The Au ho s. Published by Else ie B.V. This is an open access a icle unde he CC BY-NC-ND license (h p://c ea i ecommons.o g/licenses/by- nc-nd/4.0/). Base-age in a ian models o p edic ing indi idual ee accumula ed annual esin yield using wo apping me hods in ma i ime pine (Pinus pinas e Ai .) o es s in no h-wes e n Spain ´ Osca L´ opez-´ Al a ez * , Luis F anco-V´ azquez, Manuel Ma ey-Pe ez Resea ch G oup PROePLA GI-1716, Depa men o Plan P oduc ion and Enginee ing P ojec s, Highe Poly echnic School o Enginee ing, Campus Te a, Uni e si y o San iago de Compos ela, 27002 Lugo, Spain ARTICLE INFO Keywo ds: Pine esin P edic ing Algeb aic di e ence app oach Dynamic equa ion Be alan y-Richa ds ABSTRACT In sou he n Eu ope, especially in Spain and Po ugal, ma i ime pine esin is one o he main non- imbe o es p oduc s. A e su e ing a c isis a he end o he 20 h cen u y, i is cu en ly a g owing sec o . In Spain, depending on he a ea, he managemen o pine o es s is one o he pilla s o he na ional bioeconomy. In addi ion o imbe p oduc ion, hese o es s may be o ien ed owa ds esin p oduc ion only, o esin p oduc ion as a complemen a y ac i i y o imbe p oduc ion. In bo h cases, as in any sec o , i is essen ial o ha e ools o manage and an icipa e p oduc ion, especially in he new con ex o he bioeconomy. Fo his eason, he aim o his s udy is o de elop a dynamic model o es ima e he accumula ed esin yield du ing he esin p oduc ion season. Fo his s udy, 180 ees om h ee plo s loca ed in he no hwes o he Ibe ian Peninsula we e esin apped using wo ex ac ion me hods (non-mechanized and mechanized ci cula no ching) and s imulan pas es. Fou base models we e used om which eigh equa ions we e de i ed using ADA and GADA echniques. The mos e icien equa ions, bo h o modelling wi h he ain da a and o p edic ion wi h he es da a, we e hose de i ed om he Be alan y-Richa ds model. The RRMSE was 23% o he non-mechanised me hod and 29% o he mechanised ci cula me hod. The esul s o his s udy make i possible o add he cumula i e annual esin yield o ma i ime pine o he p ocesses ha he Be alan y-Richa ds equa ion is capable o modelling. Fu he mo e, he g ea e sa ili y o hese models will be o g ea use o he o es manage in op imising he annual ha es ing season as well as o he scien i ic communi y. 1. In oduc ion Th oughou his o y, he esou ces gene a ed by o es s ha e p o- ided di e en ypes o goods ele an o socie y (Sheppa d e al., 2020). The social and en i onmen al bene i s associa ed wi h mul i- unc ional o es managemen p ac ices ha e an impo an impac , especially in u al a eas (Lo i´ c e al., 2022). Wo ldwide, he main o es esou ce in e ms o olume and e enue is imbe , bu o he commodi ies such as non- imbe o es p oduc s (NTFPs) a e also aluable (Sa deshpande and Shackle on, 2019). Acco ding o FAO (2015), he NTFPs a e “goods de i ed om o es s ha a e angible and physical objec s o biological o igin o he han wood” like ood addi i es, ib es, esins o gums. The high added alue o NTFPs makes hem a pe ec complemen o imbe p oduc ion, p o iding an ex a sou ce o income be ween imbe ha - es s (He n´ andez-Rod íguez e al., 2017). One such bene i , de i ed om hei na u e as a biop oduc , is hei ole in mi iga ing he ca bon oo p in and educing he use o ossil-based ma e ials (Solomon, 2016). This is due o he abili y o some o hem ( ib es, gum, o esin) o ix CO 2 and o eplace pe oleum de i a i es (Demko and Macha a, 2022). In Eu ope, acco ding o Vacik e al. (2020), he economic impo ance o plan s-de i ed NTFPs emains modes , wi h he epo ed alue o NTFPs eaching 1.6 billion eu os. Da a o s a is ics a e di icul o ob ain, since a la ge p opo ion o NTFPs a e in ended o sel - consump ion (Lo i´ c e al., 2020; Winkel e al., 2022). Despi e his, he e a e Eu opean coun ies ha dedica e exclusi e esou ces o his ype o ac i i ies, such as F ance, Po ugal, and Spain, which a e commi ed o NTFPs wi h he highes added alue (Wol slehne e al., 2019). In Spain, he alue o NTFPs in ela ion o imbe a ies be ween 10 % and 50 %, depending on he sou ce (Díaz-Bal ei o e al., 2020). * Co esponding au ho . E-mail add esses: [email p o ec ed] (´ O. L´ opez-´ Al a ez), [email p o ec ed] (L. F anco-V´ azquez), [email p o ec ed] (M. Ma ey-Pe ez). Con en s lis s a ailable a ScienceDi ec Fo es Ecology and Managemen jou nal homepage: www.else ie .com/loca e/ o eco h ps://doi.o g/10.1016/j. o eco.2023.121501 Recei ed 19 June 2023; Recei ed in e ised o m 7 Sep embe 2023; Accep ed 13 Oc obe 2023 Fo es Ecology and Managemen 549 (2023) 121501 2 Acco ding o MITECO (2022), he main Spanish NTFPs p oduc s a e u les and edible ungi (106.6 million eu os), co k (83.6 million eu os), ches nu (18 million eu os), pine esin (10.3 million eu os) and pine nu s (2.6 million eu os). Pine esin is he main de ensi e ba ie agains ex e nal a acks caused by pes s and diseases o pine ees (Rissanen e al., 2019). Spanish pine esin p oduc ion has expe ienced ups and downs h oughou his o y. Since he 1980 s, esin apping ac i i y has los impo ance un il i p ac ically ceased a he beginning o his cen- u y (Soli˜ no e al., 2018). Bu since 2000 s i has been g owing in impo ance again wi hin he na ional o es y sec o and is expec ed o con inue o g ow in he coming yea s (G´ omez-Ga cía e al., 2022). This is e lec ed in he inc ease in publica ions on he subjec in some egions o Spain in ecen yea s, as in he case o Galicia (Ga cía-M´ eijome e al., 2023; L´ opez-´ Al a ez e al., 2023; Touza e al., 2021; V´ azquez-Gonz´ alez e al., 2021; Zas e al., 2020a,b). This g ow h is due o i s mul iple in- dus ial uses, he en i onmen al bene i s associa ed wi h ex ac ion and he boos he sec o is ecei ing om public and p i a e ini ia i es (Soli˜ no e al., 2018; Oono e al., 2020). The p og essi e decline o he sec o in he pas cen u y has led o a lack o echnical expe ise and knowledge gene a ion, which does no co espond o o he sec o s such as he wood sec o (G´ eno a e al., 2014; L´ opez-´ Al a ez e al., 2023a). In addi ion, exis ing da abases ha e p oblems in ob aining comple e da a se ies (Sainz e al, 2010; Calama e al., 2020). The complexi y o collec ing hese da a se ies is due o he la ge numbe o in e media e measu emen s ha need o be made and he ac ha i is a e y labo ious p ocess. In o de o make he sec o mo e compe i i e, he wo k o p oduce s mus be acili a ed by p o iding hem wi h di e en ypes o esou ces. Rega ding issues ela ed o he mechaniza ion o apping p ocess, a ious au ho s ha e compa ed he pe o mance be ween di e en mechanized and non- mechanized ex ac ion me hods (Rod íguez-Ga cía e al., 2016; L´ opez- Al a ez e al., 2023b; Ga cía-M´ eijome e al. 2023) and be ween di e en s imulan pas es used (Neis e al., 2018; Micha ila e al., 2021). In e ms o planning ools, he e is a lack o s a is ically obus models ha can p edic inal p oduc ion be o e o once he apping season has s a ed. This ype o modelling is essen ial o mul i unc ional o es esou ce planning, as i allows he s ock and expec ed alue o esinous pine o be es ima ed. This lack o obus models may be due o he complexi y o iden i- ying he explana o y a iables o o al p oduc ion, as i depends on he in e ac ion o mul iple in e ela ed a iables (Zas e al., 2020a). The pape by L´ opez-Al a ez e al. (2023b) shows he impo an biome ic and p oduc ion a iabili y o pines wi hin and be ween plo s, and p e- ious s udies mainly ela e esin yield in empe a e zones o wo en i- onmen al ac o s, empe a u e and wa e s ess (Zas e al., 2020b). These esul s a e in line wi h he ma kedly seasonal beha iou o annual esin p oduc ion, wi h he highes p oduc ions ob ained in he wa me and less ainy mon hs (Rod íguez-Ga cía e al., 2015). This pa e n causes he indi idual ee accumula ed esin p oduc ion cu e in he season o ollow a sigmoidal cu e, a pa e n widely s udied in mul iple ields (Cagla e al., 2018). Fo example, ee and s and g ow h beha es simila ly when ep esen ed o e ime (Pio esan and Biondi, 2021). Fo es e s use he cons uc ion o si e index cu es o model his end o e ime. They use h ee main echniques: guide cu e, pa ame e p edic ion and algeb aic di e ence app oach (Clu e e al., 1983). The las me hodology has he ad an age ha i allows he de elopmen o dynamic equa ions de i ed om de eloping age-in a ian models om epea ed measu emen s aken o e ime (Jo dan e al., 2006). The main me hods o de eloping hese equa ions a e he algeb aic di e ence app oach (ADA), o malized by Bailey and Clu e (1974), and he gene alized algeb aic di e ence app oach (GADA), in oduced by Cieszewski and Bailey (2000). The ADA me hod has he limi a ion ha he models de i ed om i only allow he c ea ion o amilies o anamo phic cu es o hose wi h a single asymp o e (Cieszewski, 2001). To sol e he limi a ions o he ADA me hod, he GADA me hodology was de eloped, which allows he c ea ion o amilies o polymo phic cu es wi h mul iple asymp o es (Cieszewski, 2002). In GADA, o ob ain hese amilies o cu es, he base equa ions a e expanded acco ding o se e al cha ac e is ics ha a e oo ed on di e en g ow h heo ies (Di´ eguez- A anda e al., 2005). Acco ding o hese heo ies, he cu es should s a om he o igin, be polymo phic, ha e mul iple asymp o es, and ha e he same es ima ed alue be ween p edic ion and e e ence age (Cies- zewski and Bailey, 2000). The aim o his pape is o e alua e he easibili y o modelling he indi idual ee annual accumula ed esin p oduc ion de eloping an in a ian model o pine o es s in no hwes e n Spain. To ca y ou he main objec i e, he pape has been s uc u ed in o wo pa s: (1) e al- ua ion o he alidi y o he ADA and GADA me hods o modelling he indi idual ee accumula ed annual esin yield and (2) he o ecas ing pe o mance o models. 2. Ma e ial and me hods 2.1. S udy a ea and da a The da a used o ca y ou he s udy we e limi ed due o he di icul y o collec ion and he lack o p e ious s udies. As a esul , in o ma ion o he da ase came om h ee newly es ablished plo s loca ed in low- densi y pu e s ands o Pinus pinas e in he no hwes e n Spanish e- gion o Galicia (Fig. 1). The s udy a ea is cha ac e ized by mode a e annual a e age empe a u es (11–14 ◦C) and accumula ed annual ain all be ween 800 and 1500 mm. Du ing he mon hs o June o Oc obe 2021, a o al o 180 ees we e selec ed om he plo s o esin apping p ocedu es. The minimum, mean, and maximum diame e s a b eas heigh (dbh), o al heigh s (h), and esin yield (Y) alues o he ees in each plo a e p esen ed in Table 1. The a e age dbh mee s he legal equi emen s o esin p o- cedu es in he egion. As a p elimina y s ep, he dasome ic a iables we e es ed o hei lack o in luence on esin yield by calcula ing Spea man co ela ions be ween esin and dbh (co =0.18, p- alue = 0.02) and h (co =0.07, p- alue =0.37), as nei he a iable ollowed a no mal dis ibu ion acco ding o he Shapi o-Wilk es (dbh =W =0.96, p- alue =0.001; h =W =0.98, p- alue =0.03). I was also checked whe he he yields o he plo s belonged o he same popula ion. Fo his pu pose, he Shapi o-Wilk es (W =0.969, p- Fig. 1. Dis ibu ions a ea o he ma i ime pine s ands used in he s udy. ´ O. L´ opez-´ Al a ez e al. Fo es Ecology and Managemen 549 (2023) 121501 3 alue =0.001) was used o e i y ha he da a o yields as a unc ion o plo did no sa is y he assump ion o no mali y, bu he Le ene es (F = 0.246, p- alue =0.781) did sa is y he assump ion o homoscedas ici y. Ha ing e i ied his, he K uskal-Wallis es (Fig. 2) was used o es whe he he e we e s a is ically signi ican di e ences be ween he yields o each plo , which showed ha he e we e no s a is ically sig- ni ican di e ences be ween hem. Resin apping was ca ied ou using wo di e en me hods, applying each me hod o 90 ees (30 ees pe plo using each me hod, a o al o 180 ees). One was a non-mechanized a ian o he “Ame ican” ex ac ion me hod desc ibed in Rod íguez-Ga cía e al. (2014), used mainly in he Ibe ian Peninsula and F ance. The o he ex ac ion me hod was he mechanized ci cula g oo e me hod (L´ opez-´ Al a ez e al., 2023b), in which a ci cula no ch is made in he s em o he ee and a ci cula de ice is placed o collec he esin in a closed plas ic bag. The non-mechanized me hod cu a 16 cm long and 3 cm wide ec angula s ip o ba k and cambium, while he mechanized ci cula g oo e me hod used a 5 cm diame e ool and a ba e y-powe ed sc ewd i e (cu ing ci cum e ence o 15.71 cm). A new s ip was made e e y 14 days. By he end o he apping season, which las ed i e mon hs, a o al o 10 cu s had been made. Two di e en pas es we e used as a esin p oduc ion s imulan , e hephon (8 % e hephon (60 % / ), 14 % sul- phu ic acid (50 % / ), 55 % dis illed wa e , 1.7 % polyso ba e, 1 % ce yl alcohol, 4 % aseline, 5.5 % silica, 10.8 sawdus ) and ASACIF (1 % salicylic acid, 25 % sulphu ic acid (96 % / ), 5 % p opylene glycol, 19 % whea s aw, 50 % dis illed wa e ). These pas es we e applied each ime a new s ip was made. Due o incon eniences in he da a collec ion p ocess, he e was a small g oup o ees o which some in e media e yield could no be measu ed a he p ecise momen when he new cu - ing was made and we e he e o e accumula ed in he nex weighing. To i he models, pe iodic yields we e accumula ed mon hly, gi ing a o al o i e measu emen s o each ee. 2.2. Model de elopmen To model he beha iou o esin yield o e he yea , we chose o spli he da a sample by ex ac ion me hod, bu no by s imulan pas e. This da a spli ing was done because he e was a high a iabili y be ween he Table 1 Dasome ic and esin yield cha ac e iza ion o he plo s u ilized in he s udy. dbh: diame e a b eas heigh (cm); h: o al heigh (m); Y: ee esin yield (g). Plo dbh min dbh dbh max h min h h max Y min Y Y max Culle edo 21.4 31.6 45.9 12.0 16.1 21.0 407 1419.11 3359 Pan ´ on 23.1 32.6 47.9 17.0 22.5 27.0 642 1564.47 2862 Godos 25.5 39.5 54.4 17.3 20.4 24.1 339 1492.27 2974 Fig. 2. S a is ical di e ences es be ween he esin yields o each plo . The s a is ical es s be ween g oups we e pe o m wi h s a is ical so wa e R (R Co e Team, 2022) and he “ggs a splo ” package (Pa il, 2021). ´ O. L´ opez-´ Al a ez e al. Fo es Ecology and Managemen 549 (2023) 121501 4 yields in he da a se and s a is ically signi ican di e ences be ween he ex ac ion me hods, as can be seen in Fig. 3. The e a e no s a is ical di e ences be ween he s imulan pas es (Fig. 3). 2.2.1. Models es ed Se e al models ha e been es ed, including bo h ADA and GADA models de i ed om di e en base equa ions. The e o e, he e a e equa ions wi h one o mo e a iable pa ame e s depending on he speci ic si e. The gene al no a ion (Di´ eguez-A anda e al., 2005; S´ anchez-Gonz´ alez e al., 2008; Panik, 2014) used o de elop he models has been, a 1 , a 2 ,…, a n o deno e he pa ame e s in he base models, and b 1 , b 2 ,…, b n o he pa ame e s ha we e i ed in he di e en o mu- la ions o he base equa ions. Hence, he gene al o mula ion o he de i ed models has he o m Y = ( 0 , 1 , Y 0 , b 1 , b 2 , …, b n ). These ypes o equa ions ha e mul iple uses, such as desc ibing he g ow h o ees and hei pheno ypic elemen s, o he g ow h o co k hickness, as men ioned be o e. Al hough he e a e o he uses besides hose pu ely o es y, like i s use in ishe ies (Flinn and Midway, 2021), mussels (Fuen es-San os e al., 2017) o chicken (Ma a-Es ada e al., 2020) g ow h models. Di e en ADA and GADA o mula ions om base equa ions we e es ed, bu inally he selec ed base equa ions ha achie ed con e gence (Table 2) we e he Be alan y-Richa ds (Be alan y, 1949, 1957; Richa ds, 1959), Ko (ci ed in Lundq is (1957)), Hoss eld (Hoss eld, 1822) and Weibull (Weibull, 1951; Yang e al., 1978). A o al o eigh equa ions de i ed om he p e ious base equa ions we e i ed. Fou o hem (Eqs. E1 o E4) ha e been based on he Be alan y-Richa ds model. O which, h ee (Eqs. E1 o E3) we e de i ed by he ADA me hodology, hus possessing only one si e-speci ic pa ame e . Equa ion E1 gene a es amilies o anamo phic cu es wi h mul iple asymp o es, whe eas Equa ions E2 and E3 gene a e polymo phic cu es wi h a single asymp o e (Manso e al., 2021). The esul o de i ing he Be alan y- Richa ds equa ion wi h he GADA me hodology is Eq. E4, i has mo e han one pa ame e ha is si e-speci ic and allows gene a ing amilies o polymo phic cu es wi h mul iple asymp o es (P ada e al., 2019). In he case o Eqs E5 o E8, which a e de i ed om he a ious equa ions men ioned abo e, hese a e models ob ained using he ADA me hodol- ogy, only ha ing a si e-speci ic pa ame e . 2.2.2. Pa ame e es ima ion, model selec ion and alida ion To pe o m he i and subsequen alida ion o he models, he da a wi hou ou lie s was andomly s a i ied in o qua iles acco ding o he ime lag p edic o a iable in aining (80 %) and es (20 %). The i s se was used o adjus he model pa ame e s simul aneously (local and global pa ame e s) using he dummy a iable me hod desc ibed in Cieszewski e al. (2000), and selec he one wi h he bes goodness-o - i s a is ics. The second se was used o e alua e he pe o mance o he selec ed model in an unbiased manne . Model selec ion was pe o med by nume ical and g aphical analysis o he esiduals o he i s. To e alua e he pe o mance o he models on he aining da a, goodness-o - i s a is ics such as. pseudo-R 2 (R 2 ), calcula ed acco ding o Schabenbe ge and Pie ce (2001) using he esidual sum o squa es and he o al co ec ed sum o squa es, Roo Mean Squa e E o (RMSE) and Akaike’s in o ma ion c i e ion (AIC) we e used. The p edic i e pe o mance o he models was e alua ed using he Rela i e Roo Mean Squa e E o (RRMSE). To p o ide addi ional model de ails, we ha e included he bias ( he a e age o he p edic ion esiduals) and MAE ( he a e age o he absolu e alues o he p edic ion esiduals). The s a is ical so wa e R (R Co e Team, 2022) was used o he calcula ions. The “ sample” package (F ick e al., 2022) was used o andomly sepa a e he sample in o es and aining se s. The “s a s” package (R Co e Team, 2022) was used o pe o m he non-linea models i ing. The ollowing unc ions and packages we e used o calcula e he goodness-o - i s a is ics and measu e he accu acy o he models on he es da a: he “aomisc” package (Ono i, 2020) o R 2 , he “Me ics” package (Hamne and F asco, 2018) o RMSE, he “s a s” package (R Co e Team, 2022) o AIC, and he “me ica” package (Co endo e al., 2022) o RRMSE. 3. Resul s 3.1. Accumula ed esin yield pa e n The smoo hed end o he esin yield ob ained a indi idual ha es s du ing he collec ion o esin apping da a o each o he ex ac ion me hods, ollows he pa e ns shown in Fig. 4. In he non-mechanized Fig. 3. S a is ical di e ences be ween a) he ex ac ion me hods and b) and c) he s imulan pas es acco ding o he ex ac ion me hod. The s a is ical es s be ween g oups we e pe o m wi h s a is ical so wa e R (R Co e Team, 2022) and he “ggs a splo ” package (Pa il, 2021). ´ O. L´ opez-´ Al a ez e al. Fo es Ecology and Managemen 549 (2023) 121501 5 apping me hod, which was he mos p oduc i e in e ms o o al esin yield, he immedia e esin yield be ween s ips was inc eased om he beginning o he season eaching he highes alues in he la e summe mon hs and hen dec eased un il he end o he apping season. In con as , he mechanised ci cula g oo e me hod showed a mo e consis en pa e n. I s end inc eased om he s a o he esin ex ac ion campaign un il 30 days a e he s a , emained s able un il 40 days be o e he end, and hen dec eased un il he end o he season. Consequen ly, he ci cula me hod did no show as p onounced a peak as he non-mechanised me hod. As a consequence o he phenomenon shown in Fig. 4, when he cumula i e yields o e he esin apping season o indi idual ees a e plo ed as a unc ion o esin ex ac ion me hod, hey ollow a sigmoidal pa e n (Fig. 5). I can also be seen ha he slope o he cu es o he non- mechanised me hod is g ea e han he slope o he mechanised ci cula me hod, which is a consequence o he same phenomenon. 3.2. Final models The es ima ed pa ame e s, s anda d e o s o he es ima es, goodness-o - i s a is ics and pe o mance on alida ion da a a e shown in Table 3 and Table 4. All he i ed pa ame e s we e s a is ically sig- ni ican a a 1 % le el. 3.3. Non-mechanized apping me hod All he i ed models exhibi ed sa is ac o y pe o mance a e adjus ing o cumula i e esin yield (Table 3). Upon e alua ion using he aining da a, he R 2 alues ob ained anged om 0.818 o 0.873, depending on he speci ic model. The RMSE alues a ied om 246.980 o 300.758, while he AIC alues anged om 17964.02 o 18474.63. As a measu e o he good pe o mance o he models, we looked o hose wi h he lowes RMSE and AIC alues. On he es da a, he models demons a ed ai RRMSE alues below 30 % (Li e al., 2013; Despo o ic e al., 2016), wi h none exceeding 27.018 % and he lowes being 23.189 %. The bias ange was be ween −24.34 and 19.74 g while he MAE was be ween 164.38 and 177.78 g. Du ing he aining phase, he op pe o me s in e ms o R 2 , RMSE and AIC we e h ee models de i ed om he Be alan y-Richa ds model, ollowed by one model de i ed om he Ko , Hoss eld, and Weibull equa ions. Con e sely, he leas e ec i e models we e he anamo phic and mul iple asymp o ic models de i ed om he Be alan y-Richa ds equa ion, as well as he Ko model when sol ing o he “a 1 ″ pa ame e . Simila esul s we e ob ained when assessing he p edic i e pe o mance o he models on he es da a. Among he models, equa ion E3, which co esponds o he ADA model when sol ing o he “a 3 ″ pa ame e in he Be alan y-Richa ds equa ion, exhibi ed he highes R 2 alue, as well as he lowes RMSE and AIC alues. When assessing he p edic i e pe o mance o his model on he es da a, an RRMSE o 23 %, a bias o −19.48 g and a MAE o 170.38 g we e ob ained. Fu he mo e, g aphical analysis o he e- siduals o he E3 model e ealed no e idence o he e oscedas ici y o au oco ela ion be ween he esiduals (Fig. 6). To ep esen he e olu ion o he cumula i e p oduc ions gene a ed by he model (Fig. 6), a ime ( ) o 56 days om he beginning o he apping campaign was selec ed. This pe iod was chosen as a e e ence since i is wo mon hs a e he s a o he campaign and one mon h be o e he hal way poin o he campaign, allowing o an ea ly es i- ma ion o he inpu s o he season. To ep esen he p oduc ion cu es, he di e ence be ween he maximum and minimum cumula i e p o- duc ion in he ime was di ided by ou , and his amoun was added Table 2 Base models and ADA and GADA o mula ions conside ed. Base model Pa ame e ela ed o si e Solu ion o X wi h ini ial alues ( 0 , Y 0 ) Dynamic equa ion Be alan y-Richa ds Y=a1⋅(1−e−a2⋅ )a3 a 1 =X X0=Y0 (1−e−b2⋅ 0)b3 Y=Y0⋅(1−e−b2⋅ 1−e−b2⋅ 0)b3 (E1) a 2 =X X0= −ln(1− (Y0/b1)1/b3) 0 Y=b1⋅⎛ ⎝1−(1−(Y0 b1)1/b3) / 0⎞ ⎠ b3 (E2) a 3 =X X0=ln(Y0/b1) ln(a−e−b2⋅ 0)Y=b1⋅(Y0 b1) ln(1−e−b2⋅ ) ln(1−e−b2⋅ 0) (E3) a 1 =e x a 3 =b 2 +b 3 /X X0=1 2⋅(lnY0−b2⋅L0+ (lnY0−b2⋅L0)2−4⋅b3⋅L0 √) wi h L0=ln(a−e−b1⋅ 0) Y=Y0⋅(1−e−b1⋅ 1−e−b1⋅ 0)b2+b3/X0 (E4) Ko Y=a1⋅e−a2⋅ −a3 a 1 =X X0=Y0 e−b2⋅ −b3 0 Y=Y0⋅e−b2⋅ −b3 e−b2⋅ −b3 0 (E5) a 2 =X X0= − ln(Y0/b1)/ −b3 0 Y=a1⋅(Y0 b1)( 0 )b3 (E6) Weibull Y= a1⋅(1−e(− a2⋅ a3)) a 2 =X X0= − ln(1− (Y0/b1) )/ b3 1 Y=b1⋅(1− (1−Y0/b1))( / 0)b3 (E7) Hoss eld Y=a1 1+a2⋅ −a3 A 2 =X X0= −b3 0⋅(b1/Y0−1)Y=b1 1−(1−b1 Y0)⋅( 0 )b3 (E8) Fig. 4. Resin yields o indi idual g oo es du ing he apping season o he wo me hods, non-mechanized apping me hod (con inuous line) and mechanized ci cula g oo e me hod (dashed line). ´ O. L´ opez-´ Al a ez e al. Fo es Ecology and Managemen 549 (2023) 121501 6 om he minimum o he maximum p oduc ion in he pe iod . The ep esen ed cu es a e hose o 187.00, 410.25, 633.50, 856.75, and 1080.00 accumula ed g ams a ime . Fig. 8 plo s he RRMSE ob ained by he model when making p edic ions on es da a as a unc ion o ime. The model achie es an RRMSE o less han 30 % o o ecas s up o 90 days (Fig. 8), which is e y good o o ecas s o less han one mon h, good o o ecas s o wo mon hs and ai o o ecas s o h ee mon hs acco ding o Li e al. (2013). This s a is ic ises abo e 30 % o o ecas s longe han 90 days. 3.4. Mechanized ci cula g oo e me hod The i s o he da a ob ained wi h he mechanized ci cula g oo e me hod (Table 4), show a ou able esul s bu mo e modes compa ed o he non-mechanized apping me hod. The R 2 alues achie ed wi h he aining da a ange om 0.695 o 0.82. The RMSE alues and AIC alues ange om 212.512 o 276.569 and om 16844.4 o 17496.84, espec i ely. When assessing he models on es da a, he RRMSE was again ob ained alues below 30 %. The bias a ied om 0.98 o 41.45 g, and he MAE anged om 150.95 and 199.56 g. Consis en ly, equa ions de i ed om he Be alan y-Richa ds model demons a e supe io goodness-o - i s a is ics, yielding he mos op imal i o he da a. Con e sely, he poo es i s a e obse ed wi h he same equa ions as in he non-mechanized apping me hod case, speci - ically he anamo phic and mul iple asymp o ic o m o he Be alan y- Richa ds equa ion, as well as he Ko model by sol ing he “a 1 ″ pa ame e . The E4 o GADA o m o he Be alan y-Richa ds model is ound o Fig. 5. Accumula ed esin yield du ing esin apping season o bo h esin apping me hods, non-mechanized apping me hod (g aphic a he op) and mechanized ci cula g oo e me hod (g aphic a he bo om). ´ O. L´ opez-´ Al a ez e al. Fo es Ecology and Managemen 549 (2023) 121501 7 be he equa ion ha achie es he highes le el o s a is ical goodness o i wi h he aining da a (Table 4). G aphical analysis o he esiduals o he E4 model indica es he absence o any obse able he e o- scedas ici y o au oco ela ion among he esiduals (Fig. 7). When es ed agains he alida ion da a, he equa ion displayed he lowes RRMSE alue o 29 %, a bias o 0.98 and a MAE o 150.95 g. To ep esen he cu es o he adjus men s made wi h he equa ion on he o iginal da a, i has been done in he same way as in he p e ious Table 3 Pa ame e es ima es and goodness-o - i s a is ics o non-mechanized apping me hod. Model T ain Tes Pa ame e Es ima e S anda d e o R 2 RMSE AIC RRMSE (%) Bias MAE E1 b 2 −0.003 0.001 0.818 295.901 18432.44 26.370 16.09 177.78 b 3 0.984 0.037 E2 b 1 4180.555 16.090 0.856 263.305 18129.92 24.806 −14.66 167.41 b 3 1.535 0.028 E3 b 1 15767.82 2.216x103 0.873 246.980 17964.02 23.189 −19.48 170.38 b 2 1.497x10−3 2.148x10−4 E4 b 1 1.608x10−3 1.425x10−4 0.873 247.330 17969.69 23.415 −24.34 170.85 b 2 −1.029x104 1.284x103 b 3 9.908x104 1.228x104 E5 b 2 17 3.582 0.812 300.758 18474.63 27.018 19.74 175.57 b 3 0.096 0.034 E6 b 1 3.447x104 7.237x103 0.857 261.883 18115.89 24.198 −11.32 164.38 b 3 0.318 0.019 E7 b 1 3.768x103 1.286x102 0.854 264.460 18141.26 25.047 −14.12 167.66 b 3 1.398 1.631x10−2 E8 b 1 5.233x103 2.449x102 0.855 263.686 18133.66 24.827 −14.19 167.72 b 3 1.465 2.202x10−2 Table 4 Pa ame e es ima es and goodness-o - i s a is ics o mechanized ci cula g oo e me hod. Model T ain Tes Pa ame e Es ima e S anda d e o R 2 RMSE AIC RRMSE (%) Bias MAE E1 b 2 0.005 0.001 0.696 276.414 17495.44 41.971 41.45 199.43 b 3 1.247 0.059 E2 b 1 2.845x103 64.180 0.770 236.704 17110.2 33.698 12.67 170.69 b 3 1.589 0.029 E3 b 1 4.935x103 3.447x102 0.811 217.614 16901.33 29.844 4.29 156.77 b 2 3.603x10−3 3.076x10−4 E4 b 1 1.077x10−2 9.118x10−4 0.820 212.512 16844.4 29.244 0.98 150.95 b 2 −13.780 1.543 b 3 1.197x102 12.040 E5 b 2 11.931 0.629 0.695 276.569 17496.84 41.916 44.63 199.56 b 3 0.189 0.046 E6 b 1 5.885x103 4.574x102 0.817 214.433 16864.75 29.437 0.83 152.23 b 3 0.541 2.092x10−2 E7 b 1 2820.459 67.679 0.763 244.104 17186.67 35.540 21.69 178.76 b 3 1.357 0.017 E8 b 1 3.120x103 92.570 0.781 234.208 17083.87 33.020 12.54 167.14 b 3 1.563 0.024 Fig. 6. The le igu e was he age-dependen dynamic model E3 o non-mechanized no ching and he accumula ed esin yield cu es o 187.00, 410.25, 633.50, 856.75, and 1080.00 g a he ime o 56 days. The igh igu e is he Residuals s Residuals i +1 o non-mechanized no ching model E3. ´ O. L´ opez-´ Al a ez e al. Fo es Ecology and Managemen 549 (2023) 121501 8 case. In his case, he cu es ep esen ed a ime a e hose o 55.00, 333.75, 612.50, 891.25, and 1170.00 g (Fig. 7). The model ob ained an RRMSE lowe han 35 % when he p edic ion lag was 85 days (Fig. 8), also educing he mean p edic ion e o a he end o he apping season below ha ob ained by he model i ed o he o he ex ac ion me hod. In his case, o ecas s o less han one mon h could be desc ibed as good and hose o less han wo and a hal mon hs as ai (Li e al., 2013; Despo o ic e al., 2016). 4. Discussion This s udy ep esen s he i s case o modelling indi idual ee accumula i e esin yield h oughou he esin p oduc ion season. I demons a es ha me hods simila o hose used o model di e en na u al g ow h p ocesses can be success ully applied. The unp ece- den ed na u e o his p oduc ion modelling poses a challenge o he discussion, gi en he lack o p e ious s udies on esin o es ablish a compa a i e amewo k. Fi ing models desc ibing he end o cumu- la i e p oduc ion du ing he NTFP ha es season is challenging due o he high a iabili y in indi idual ee p oduc ion (G´ omez-Ga cía e al., 2022; L´ opez-Al a ez e al., 2023b) and he in luence ha en i onmen changes ha e on he yield o his enewable p oduc (V´ azquez-Gonz´ alez e al., 2021b). Rega dless o he ex ac ion me hod, he equa ions de i ed om he Be alan y-Richa ds model wi h ADA and GADA me hodologies we e he ones ha p esen ed he bes pe o mance. Likewise, when compa ing he s a is ics o he models based on he ex ac ion me hods, he model o he non-mechanized me hod ob ained be e goodness-o - i han he mechanized me hod. The ob ained esul s sugges ha non-mechanized apping may lead o a mo e p onounced yield inc ease in indi idual g oo es owa ds he middle and end o he season when compa ed o he mechanized ci cula g oo e me hod (Fig. 4). This beha io , de i ed mainly om he apping me hod (L´ opez-´ Al a ez e al., 2023b), accen ua es he p esence o an in lec ion poin in he cumula i e p oduc ion cu e, making i mo e simila o g ow h p ocesses. In he non-mechanized apping me hod, he equa ions ha ob ained he highes R 2 in he aining phase we e hose de i ed om he Be alan y-Richa ds equa ion by ADA and GADA me hodology (E3 and E4, espec i ely). The simila i y in desc ibing he end o he da a be ween he ADA and GADA o ms o he Be alan y- Richa ds model is common in o he o es species g ow h modelling wo ks (T im e al., 2020; Wang e al., 2020). The GADA o m gene ally ob ains be e goodness-o - i s a is ics, as i can gene a e polymo phic cu es wi h mul iple asymp o es. This s a emen is e i ied in he case o he non-mechanized ex ac ion me hod, as he GADA o mula ion ob- ained he bes s a is ics, simila o he epo ed by S´ anchez-Gonz´ alez Fig. 7. The le igu e was he age-dependen dynamic model E4 o mechanized ci cula g oo e and he accumula ed esin yield cu es o 55.00, 333.75, 612.50, 891.25, and 1170.00 g a he ime o 56 days. The igh igu e is he Residuals s Residuals i +1 o mechanized ci cula g oo e model E4. Fig. 8. Rela i e oo mean squa e e o (RRMSE) by lag o p edic ion o mechanized ci cula g oo e (ci cles) and non-mechanized no ching ( iangles). ´ O. L´ opez-´ Al a ez e al. Fo es Ecology and Managemen 549 (2023) 121501 9 e al. (2008) when modeling co k g ow h. The goodness-o - i s a is ics ob ained in bo h i s a e simila o hose ob ained when using he GADA o m o he Be alan y-Richa ds equa ion o model p oduc ion o e ime o long o a ion species such as he peduncula e oak (G´ omez-Ga cía e al., 2015), and sligh ly lowe han hose ob ained when modeling sho o a ion species (Di´ eguez-A anda e al., 2005; Seki and Sakici, 2017) o co k p oduc ion (S´ anchez-Gonz´ alez e al., 2008). The good pe o mance o he models in e ms o goodness-o - i , as ob ained in his g ound-b eaking s udy, e i ies ha hese echniques can be used o model esin p oduc ion du ing he season. To make hese models a alid ool o modelling cumula i e esin p oduc ion, bo h in he ield o ac i i y o o es wo ke s, in he esin apping indus y and in he esin p ocessing indus y, i would be necessa y o de elop his p elimina y s udy u he and o ca y ou a la ge sample, ep esen a- i e o he whole en i onmen al g adien and including all he possible casuis y unde which esin pine o es s can be ound. The e o e, based on he abo e, desc ibing he indi idual ee accumula ed esin p oduc- ion o P. pinas e can be included among he uses o he on Be alan y model (Be alan y, 1949, 1957). The p edic ions made wi h he es da a yielded a e age RRMSE alues below 30 %, ega dless o he ex ac ion me hod. These a e age alues can be classi ied as ai acco ding o Despo o ic e al. (2016) and Li e al. (2013). Simpli ying and no conside ing he pa icula lag o each p edic ion, he models ha e an a e age accu acy o 76.81 % and 70.76 %, o he non-mechanized and mechanized me hods, espec- i ely. These accu acy pe cen ages a e sligh ly lowe han hose ach- ie ed by he g ow h models used in dend ome y (P ada e al., 2019) and a e mo e simila o hose ob ained by S´ anchez-Gonz´ alez e al. (2008) modeling co k g ow h. When e alua ing he accu acy pe cen age as a unc ion o he lag o he p edic ions, i was ound ha a p edic ion made wi h an 80-day lag, which is usually hal o he esin- apping campaign, had an e o a e o app oxima ely 35 %. These accu acy pe cen ages as a unc ion o he lag a e sligh ly lowe han hose ach- ie ed by dend ome ic models (Gea-Izquie do e al., 2008). The lowe p ecision alues may be due o he small sample size o he high a i- abili y o esin p oduc ion o e ime due o dasome ic and en i on- men al a iables, as da a om dend ome ic a iables end o show less a iabili y in p oduc ion be ween measu emen pe iods. To acili a e scien i ic esea ch, one po en ial unc ionali y o he esul s ob ained is hei use ulness o eco e missing da a om he empo al weigh se ies. Collec ing in e media e weigh da a is a di icul ask (Calama e al., 2011), and usually he e a e da a ha we e no measu ed a he p ecise momen when each new s ip was done bu we e accumula ed and weighed a e wo o mo e apping pe iods. In his way, he accumula ed p oduc ion is accoun ed o , bu i is no known which p oduc ion belongs o each pe iod. By using his ype o models, he p oduc ion cu es can be ec ea ed, and he missing in e media e da a can be eco e ed. They a e pa icula ly ele an o esea che s who wish o gene a e and moni o p oduc ion cu es and ha e missing da a in he p oduc ion se ies. The good p edic i e capaci y o he models makes hem a alid ool o he mul i unc ional managemen o pine o es s, which is easible as long as he owne s a e p o ided wi h he app op ia e ools o ca y i ou (Sabas ian e al., 2019). By using hese models, o es manage s will be able o p edic inpu s be o e he esin ha es ing season ends, allowing hem o op imize he du a ion o he campaign o maximise he annual income om esin sales. This makes esin ex ac ion mo e a ac i e o o es p oduce s, c ea ing g een jobs and p omo ing he bioeconomy by encou aging he p oduc ion o a biop oduc such as pine esin. 5. Conclusion Wi h he inc easing impo ance o esin p oduc ion in he NTFPs and o es y sec o , he e is a need o u he de elop and expand ou knowledge. This wo k demons a es o he i s ime ha age-in a ian models a e a alid ool o he managemen o pine o es s unde esin apping p ocedu es. Acco ding o R 2 , RMSE, AIC, RRMSE, Bias and MAE, he ADA and GADA o mula ions o he Be alan y-Richa ds model pe o med bes in modelling cumula i e esin p oduc ion. This implies he inclusion o accumula i e esin yield o indi idual ees as one o he many phe- nomena ha can be modelled wi h hese equa ions. Fu he mo e, he abili y o accu a ely p edic inal p oduc ion once he campaign has s a ed using hese models ep esen s a signi ican ad ance in mul i- unc ional o es managemen . In addi ion, his ype o model applied o cumula i e annual esin p oduc ion can p o ide a basis o econ- s uc ing missing da a se ies wi h a high deg ee o accu acy. To u he ad ance his sec o , public o p i a e ini ia i es should p o ide suppo and ollow he example o o es g ow h models by es ablishing a ne wo k o pe manen plo s o pa ame ising na ional-scale esin p o- duc ion a iables and de eloping mo e obus models. Enabling ad ancemen in he e ec i e join managemen o wood esou ces and esin p oduc ion, acili a ing he c ea ion o echnical–economic models and he implemen a ion o inno a i e, cu ing-edge me hodologies in cons uc ing said models. Founding and acknowledgemen s This wo k was suppo ed by he Spanish Go e nmen (“ACREMA”, MAPA/AEI-Ag i/FEADER, UE) [O00000226e2000043659], he Gali- cian Go e nmen (Xun a de Galicia) wi h a g an o Compe i i e Re e ence G oups ED431C-2021-27 and he p e-doc o al con ac Campus Te a-USC 2023. The au ho s hank FORESIN and CIF Lou izan, o he ield assessmen s. CRediT au ho ship con ibu ion s a emen ´ Osca L´ opez-´ Al a ez: Concep ualiza ion, Me hodology, Fo mal analysis, W i ing – o iginal d a , Visualiza ion. Luis F anco-V´ azquez: Me hodology, W i ing – e iew & edi ing. Manuel Ma ey-Pe ez: Concep ualiza ion, Me hodology, Fo mal analysis, Resou ces, W i ing – e iew & edi ing, Supe ision, P ojec adminis a ion, Funding acquisi ion. Decla a ion o Compe ing In e es The au ho s decla e ha hey ha e no known compe ing inancial in e es s o pe sonal ela ionships ha could ha e appea ed o in luence he wo k epo ed in his pape . Da a a ailabili y Pe mission o sha e he da a mus be ob ained om he ins i u ions ha unded he esea ch Re e ences Bailey, R.L., Clu e , J.L., 1974. Base-age in a ian polymo phic si e cu es. Fo . Sci. 20, 155–159. h ps://doi.o g/10.1093/ o es science/20.2.155. Be alan y, L.V., 1957. Quan i a i e laws in me abolism and g ow h. Q. Re . Biol. 32, 217–231. h ps://doi.o g/10.1086/401873. Cagla , M.U., Teu el, A.I., Wilke, C.O., 2018. Sicega : R package o sigmoidal and double-sigmoidal cu e i ing. Pee J 6, e4251. Calama R., Miina, J., de-Miguel, S., Bone , J. A., Mouni , F., Tom´ e, M., Ma ínez- Jaú egui, M., He uzo, C., Pel ola, R., Salo, K., Ku ila M., He n´ andez-Rod íguez, M., Ma ín-Pin o, P., S´ anchez-Gonz´ alez, M., 2020. Da a & models: impo ance o assessing and o ecas ing non-wood o es p oduc s in Eu ope, in: Vacik, H., Hale, M., Spiecke , H., Pe enella, D. Tom´ e, M., Non-Wood Fo es P oduc s in Eu ope. Ecology and managemen o mush ooms, ee p oduc s, unde s o y plan s and animal p oduc s. Ou comes o he COST Ac ion FP1203 on Eu opean NWFPs, BoD, No de s ed , pp. 43-78. Calama Sainz, R., Tome, M., S´ anchez-Gonz´ alez, M., Miina., J., Spanos, K., Palahi, M., 2011. Modelling Non-Wood Fo es P oduc s in Eu ope: a e iew. Fo es Sys . 3 (4), 69. ´ O. L´ opez-´ Al a ez e al.