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

Forest Ecology and Management 549 (2023) 121501 Available online 19 October 2023 0378-1127/© 2023 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). 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 ´ Oscar L´ opez-´ Alvarez * , Luis Franco-V´ azquez, Manuel Marey-Perez Research Group PROePLA GI-1716, Department of Plant Production and Engineering Projects, Higher Polytechnic School of Engineering, Campus Terra, University of Santiago de Compostela, 27002 Lugo, Spain ARTICLE INFO Keywords: Pine resin Predicting Algebraic difference approach Dynamic equation Bertalanffy-Richards ABSTRACT In southern Europe, especially in Spain and Portugal, maritime pine resin is one of the main non-timber forest products. After suffering a crisis at the end of the 20th century, it is currently a growing sector. In Spain, depending on the area, the management of pine forests is one of the pillars of the national bioeconomy. In addition to timber production, these forests may be oriented towards resin production only, or resin production as a complementary activity to timber production. In both cases, as in any sector, it is essential to have tools to manage and anticipate production, especially in the new context of the bioeconomy. For this reason, the aim of this study is to develop a dynamic model to estimate the accumulated resin yield during the resin production season. For this study, 180 trees from three plots located in the northwest of the Iberian Peninsula were resin tapped using two extraction methods (non-mechanized and mechanized circular notching) and stimulant pastes. Four base models were used from which eight equations were derived using ADA and GADA techniques. The most efficient equations, both for modelling with the train data and for prediction with the test data, were those derived from the Bertalanffy-Richards model. The RRMSE was 23% for the non-mechanised method and 29% for the mechanised circular method. The results of this study make it possible to add the cumulative annual resin yield of maritime pine to the processes that the Bertalanffy-Richards equation is capable of modelling. Furthermore, the great versatility of these models will be of great use to the forest manager in optimising the annual harvesting season as well as for the scientific community. 1. Introduction Throughout history, the resources generated by forests have provided different types of goods relevant to society (Sheppard et al., 2020). The social and environmental benefits associated with multi-functional forest management practices have an important impact, especially in rural areas (Lovri´ c et al., 2022). Worldwide, the main forest resource in terms of volume and revenue is timber, but other commodities such as non-timber forest products (NTFPs) are also valuable (Sardeshpande and Shackleton, 2019). According to FAO (2015), the NTFPs are “goods derived from forests that are tangible and physical objects of biological origin other than wood” like food additives, fibres, resins or gums. The high added value of NTFPs makes them a perfect complement to timber production, providing an extra source of income between timber harvests (Hern´ andez-Rodríguez et al., 2017). One such benefit, derived from their nature as a bioproduct, is their role in mitigating the carbon footprint and reducing the use of fossil-based materials (Solomon, 2016). This is due to the ability of some of them (fibres, gum, or resin) to fix CO 2 and to replace petroleum derivatives (Demko and Machava, 2022). In Europe, according to Vacik et al. (2020), the economic importance of plants-derived NTFPs remains modest, with the reported value of NTFPs reaching 1.6 billion euros. Data for statistics are difficult to obtain, since a large proportion of NTFPs are intended for selfconsumption (Lovri´ c et al., 2020; Winkel et al., 2022). Despite this, there are European countries that dedicate exclusive resources to this type of activities, such as France, Portugal, and Spain, which are committed to NTFPs with the highest added value (Wolfslehner et al., 2019). In Spain, the value of NTFPs in relation to timber varies between 10 % and 50 %, depending on the source (Díaz-Balteiro et al., 2020). * Corresponding author. E-mail addresses: [email protected] (´ O. L´ opez-´ Alvarez), [email protected] (L. Franco-V´ azquez), [email protected] (M. Marey-Perez). Contents lists available at ScienceDirect Forest Ecology and Management journal homepage: www.elsevier.com/locate/foreco https://doi.org/10.1016/j.foreco.2023.121501 Received 19 June 2023; Received in revised form 7 September 2023; Accepted 13 October 2023 Forest Ecology and Management 549 (2023) 121501 2 According to MITECO (2022), the main Spanish NTFPs products are truffles and edible fungi (106.6 million euros), cork (83.6 million euros), chestnut (18 million euros), pine resin (10.3 million euros) and pine nuts (2.6 million euros). Pine resin is the main defensive barrier against external attacks caused by pests and diseases of pine trees (Rissanen et al., 2019). Spanish pine resin production has experienced ups and downs throughout history. Since the 1980 s, resin tapping activity has lost importance until it practically ceased at the beginning of this century (Soli˜ no et al., 2018). But since 2000 s it has been growing in importance again within the national forestry sector and is expected to continue to grow in the coming years (G´ omez-García et al., 2022). This is reflected in the increase in publications on the subject in some regions of Spain in recent years, as in the case of Galicia (García-M´ eijome et al., 2023; L´ opez-´ Alvarez et al., 2023; Touza et al., 2021; V´ azquez-Gonz´ alez et al., 2021; Zas et al., 2020a,b). This growth is due to its multiple industrial uses, the environmental benefits associated with extraction and the boost the sector is receiving from public and private initiatives (Soli˜ no et al., 2018; Oono et al., 2020). The progressive decline of the sector in the past century has led to a lack of technical expertise and knowledge generation, which does not correspond to other sectors such as the wood sector (G´ enova et al., 2014; L´ opez-´ Alvarez et al., 2023a). In addition, existing databases have problems in obtaining complete data series (Sainz et al, 2010; Calama et al., 2020). The complexity of collecting these data series is due to the large number of intermediate measurements that need to be made and the fact that it is a very laborious process. In order to make the sector more competitive, the work of producers must be facilitated by providing them with different types of resources. Regarding issues related to the mechanization of tapping process, various authors have compared the performance between different mechanized and nonmechanized extraction methods (Rodríguez-García et al., 2016; L´ opezAlvarez et al., 2023b; García-M´ eijome et al. 2023) and between different stimulant pastes used (Neis et al., 2018; Michavila et al., 2021). In terms of planning tools, there is a lack of statistically robust models that can predict final production before or once the tapping season has started. This type of modelling is essential for multifunctional forest resource planning, as it allows the stock and expected value of resinous pine to be estimated. This lack of robust models may be due to the complexity of identifying the explanatory variables of total production, as it depends on the interaction of multiple interrelated variables (Zas et al., 2020a). The paper by L´ opez-Alvarez et al. (2023b) shows the important biometric and production variability of pines within and between plots, and previous studies mainly relate resin yield in temperate zones to two environmental factors, temperature and water stress (Zas et al., 2020b). These results are in line with the markedly seasonal behaviour of annual resin production, with the highest productions obtained in the warmer and less rainy months (Rodríguez-García et al., 2015). This pattern causes the individual tree accumulated resin production curve in the season to follow a sigmoidal curve, a pattern widely studied in multiple fields (Caglar et al., 2018). For example, tree and stand growth behaves similarly when represented over time (Piovesan and Biondi, 2021). Foresters use the construction of site index curves to model this trend over time. They use three main techniques: guide curve, parameter prediction and algebraic difference approach (Clutter et al., 1983). The last methodology has the advantage that it allows the development of dynamic equations derived from developing age-invariant models from repeated measurements taken over time (Jordan et al., 2006). The main methods for developing these equations are the algebraic difference approach (ADA), formalized by Bailey and Clutter (1974), and the generalized algebraic difference approach (GADA), introduced by Cieszewski and Bailey (2000). The ADA method has the limitation that the models derived from it only allow the creation of families of anamorphic curves or those with a single asymptote (Cieszewski, 2001). To solve the limitations of the ADA method, the GADA methodology was developed, which allows the creation of families of polymorphic curves with multiple asymptotes (Cieszewski, 2002). In GADA, to obtain these families of curves, the base equations are expanded according to several characteristics that are rooted on different growth theories (Di´ eguezAranda et al., 2005). According to these theories, the curves should start from the origin, be polymorphic, have multiple asymptotes, and have the same estimated value between prediction and reference age (Cieszewski and Bailey, 2000). The aim of this paper is to evaluate the feasibility of modelling the individual tree annual accumulated resin production developing an invariant model for pine forests in northwestern Spain. To carry out the main objective, the paper has been structured into two parts: (1) evaluation of the validity of the ADA and GADA methods for modelling the individual tree accumulated annual resin yield and (2) the forecasting performance of models. 2. Material and methods 2.1. Study area and data The data used to carry out the study were limited due to the difficulty of collection and the lack of previous studies. As a result, information for the dataset came from three newly established plots located in lowdensity pure stands of Pinus pinaster in the northwestern Spanish region of Galicia (Fig. 1). The study area is characterized by moderate annual average temperatures (11–14 ◦C) and accumulated annual rainfall between 800 and 1500 mm. During the months of June to October 2021, a total of 180 trees were selected from the plots for resin tapping procedures. The minimum, mean, and maximum diameters at breast height (dbh), total heights (h), and resin yield (Y) values for the trees in each plot are presented in Table 1. The average dbh meets the legal requirements for resin procedures in the region. As a preliminary step, the dasometric variables were tested for their lack of influence on resin yield by calculating Spearman correlations between resin and dbh (corr =0.18, p-value = 0.02) and h (corr =0.07, p-value =0.37), as neither variable followed a normal distribution according to the Shapiro-Wilk test (dbh =W =0.96, p-value =0.001; h =W =0.98, p-value =0.03). It was also checked whether the yields of the plots belonged to the same population. For this purpose, the Shapiro-Wilk test (W =0.969, pFig. 1. Distributions area of the maritime pine stands used in the study. ´ O. L´ opez-´ Alvarez et al. Forest Ecology and Management 549 (2023) 121501 3 value =0.001) was used to verify that the data of yields as a function of plot did not satisfy the assumption of normality, but the Levene test (F = 0.246, p-value =0.781) did satisfy the assumption of homoscedasticity. Having verified this, the Kruskal-Wallis test (Fig. 2) was used to test whether there were statistically significant differences between the yields of each plot, which showed that there were no statistically significant differences between them. Resin tapping was carried out using two different methods, applying each method to 90 trees (30 trees per plot using each method, a total of 180 trees). One was a non-mechanized variant of the “American” extraction method described in Rodríguez-García et al. (2014), used mainly in the Iberian Peninsula and France. The other extraction method was the mechanized circular groove method (L´ opez-´ Alvarez et al., 2023b), in which a circular notch is made in the stem of the tree and a circular device is placed to collect the resin in a closed plastic bag. The non-mechanized method cut a 16 cm long and 3 cm wide rectangular strip of bark and cambium, while the mechanized circular groove method used a 5 cm diameter tool and a battery-powered screwdriver (cutting circumference of 15.71 cm). A new strip was made every 14 days. By the end of the tapping season, which lasted five months, a total of 10 cuts had been made. Two different pastes were used as a resin production stimulant, ethephon (8 % ethephon (60 % v/v), 14 % sulphuric acid (50 % v/v), 55 % distilled water, 1.7 % polysorbate, 1 % cetyl alcohol, 4 % vaseline, 5.5 % silica, 10.8 sawdust) and ASACIF (1 % salicylic acid, 25 % sulphuric acid (96 % v/v), 5 % propylene glycol, 19 % wheat straw, 50 % distilled water). These pastes were applied each time a new strip was made. Due to inconveniences in the data collection process, there was a small group of trees for which some intermediate yield could not be measured at the precise moment when the new cutting was made and were therefore accumulated in the next weighing. To fit the models, periodic yields were accumulated monthly, giving a total of five measurements for each tree. 2.2. Model development To model the behaviour of resin yield over the year, we chose to split the data sample by extraction method, but not by stimulant paste. This data splitting was done because there was a high variability between the Table 1 Dasometric and resin yield characterization of the plots utilized in the study. dbh: diameter at breast height (cm); h: total height (m); Y: tree resin yield (g). Plot dbh min dbh dbh max h min h h max Y min Y Y max Culleredo 21.4 31.6 45.9 12.0 16.1 21.0 407 1419.11 3359 Pant´ 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. Statistical differences test between the resin yields of each plot. The statistical tests between groups were perform with statistical software R (R Core Team, 2022) and the “ggstatsplot” package (Patil, 2021). ´ O. L´ opez-´ Alvarez et al. Forest Ecology and Management 549 (2023) 121501 4 yields in the data set and statistically significant differences between the extraction methods, as can be seen in Fig. 3. There are no statistical differences between the stimulant pastes (Fig. 3). 2.2.1. Models tested Several models have been tested, including both ADA and GADA models derived from different base equations. Therefore, there are equations with one or more variable parameters depending on the specific site. The general notation (Di´ eguez-Aranda et al., 2005; S´ anchez-Gonz´ alez et al., 2008; Panik, 2014) used to develop the models has been, a 1 , a 2 ,…, a n to denote the parameters in the base models, and b 1 , b 2 ,…, b n for the parameters that were fitted in the different formulations of the base equations. Hence, the general formulation of the derived models has the form Y =f (t 0 , t 1 , Y 0 , b 1 , b 2 , …, b n ). These types of equations have multiple uses, such as describing the growth of trees and their phenotypic elements, or the growth of cork thickness, as mentioned before. Although there are other uses besides those purely forestry, like its use in fisheries (Flinn and Midway, 2021), mussels (Fuentes-Santos et al., 2017) or chicken (Mata-Estrada et al., 2020) growth models. Different ADA and GADA formulations from base equations were tested, but finally the selected base equations that achieved convergence (Table 2) were the Bertalanffy-Richards (Bertalanffy, 1949, 1957; Richards, 1959), Korf (cited in Lundqvist (1957)), Hossfeld (Hossfeld, 1822) and Weibull (Weibull, 1951; Yang et al., 1978). A total of eight equations derived from the previous base equations were fitted. Four of them (Eqs. E1 to E4) have been based on the Bertalanffy-Richards model. Of which, three (Eqs. E1 to E3) were derived by the ADA methodology, thus possessing only one site-specific parameter. Equation E1 generates families of anamorphic curves with multiple asymptotes, whereas Equations E2 and E3 generate polymorphic curves with a single asymptote (Manso et al., 2021). The result of deriving the BertalanffyRichards equation with the GADA methodology is Eq. E4, it has more than one parameter that is site-specific and allows generating families of polymorphic curves with multiple asymptotes (Prada et al., 2019). In the case of Eqs E5 to E8, which are derived from the various equations mentioned above, these are models obtained using the ADA methodology, only having a site-specific parameter. 2.2.2. Parameter estimation, model selection and validation To perform the fit and subsequent validation of the models, the data without outliers was randomly stratified into quartiles according to the time lag predictor variable in training (80 %) and test (20 %). The first set was used to adjust the model parameters simultaneously (local and global parameters) using the dummy variable method described in Cieszewski et al. (2000), and select the one with the best goodness-of-fit statistics. The second set was used to evaluate the performance of the selected model in an unbiased manner. Model selection was performed by numerical and graphical analysis of the residuals of the fits. To evaluate the performance of the models on the training data, goodness-of-fit statistics such as. pseudo-R 2 (R 2 ), calculated according to Schabenberger and Pierce (2001) using the residual sum of squares and the total corrected sum of squares, Root Mean Square Error (RMSE) and Akaike’s information criterion (AIC) were used. The predictive performance of the models was evaluated using the Relative Root Mean Square Error (RRMSE). To provide additional model details, we have included the bias (the average of the prediction residuals) and MAE (the average of the absolute values of the prediction residuals). The statistical software R (R Core Team, 2022) was used for the calculations. The “rsample” package (Frick et al., 2022) was used to randomly separate the sample into test and training sets. The “stats” package (R Core Team, 2022) was used to perform the non-linear models fitting. The following functions and packages were used to calculate the goodness-of-fit statistics and measure the accuracy of the models on the test data: the “aomisc” package (Onofri, 2020) for R 2 , the “Metrics” package (Hamner and Frasco, 2018) for RMSE, the “stats” package (R Core Team, 2022) for AIC, and the “metrica” package (Correndo et al., 2022) for RRMSE. 3. Results 3.1. Accumulated resin yield pattern The smoothed trend of the resin yield obtained at individual harvests during the collection of resin tapping data for each of the extraction methods, follows the patterns shown in Fig. 4. In the non-mechanized Fig. 3. Statistical differences between a) the extraction methods and b) and c) the stimulant pastes according to the extraction method. The statistical tests between groups were perform with statistical software R (R Core Team, 2022) and the “ggstatsplot” package (Patil, 2021). ´ O. L´ opez-´ Alvarez et al. Forest Ecology and Management 549 (2023) 121501 5 tapping method, which was the most productive in terms of total resin yield, the immediate resin yield between strips was increased from the beginning of the season reaching the highest values in the late summer months and then decreased until the end of the tapping season. In contrast, the mechanised circular groove method showed a more consistent pattern. Its trend increased from the start of the resin extraction campaign until 30 days after the start, remained stable until 40 days before the end, and then decreased until the end of the season. Consequently, the circular method did not show as pronounced a peak as the non-mechanised method. As a consequence of the phenomenon shown in Fig. 4, when the cumulative yields over the resin tapping season of individual trees are plotted as a function of resin extraction method, they follow a sigmoidal pattern (Fig. 5). It can also be seen that the slope of the curves of the nonmechanised method is greater than the slope of the mechanised circular method, which is a consequence of the same phenomenon. 3.2. Final models The estimated parameters, standard errors of the estimates, goodness-of-fit statistics and performance on validation data are shown in Table 3 and Table 4. All the fitted parameters were statistically significant at a 1 % level. 3.3. Non-mechanized tapping method All the fitted models exhibited satisfactory performance after adjusting for cumulative resin yield (Table 3). Upon evaluation using the training data, the R 2 values obtained ranged from 0.818 to 0.873, depending on the specific model. The RMSE values varied from 246.980 to 300.758, while the AIC values ranged from 17964.02 to 18474.63. As a measure of the good performance of the models, we looked for those with the lowest RMSE and AIC values. On the test data, the models demonstrated fair RRMSE values below 30 % (Li et al., 2013; Despotovic et al., 2016), with none exceeding 27.018 % and the lowest being 23.189 %. The bias range was between −24.34 and 19.74 g while the MAE was between 164.38 and 177.78 g. During the training phase, the top performers in terms of R 2 , RMSE and AIC were three models derived from the Bertalanffy-Richards model, followed by one model derived from the Korf, Hossfeld, and Weibull equations. Conversely, the least effective models were the anamorphic and multiple asymptotic models derived from the Bertalanffy-Richards equation, as well as the Korf model when solving for the “a 1 ″ parameter. Similar results were obtained when assessing the predictive performance of the models on the test data. Among the models, equation E3, which corresponds to the ADA model when solving for the “a 3 ″ parameter in the Bertalanffy-Richards equation, exhibited the highest R 2 value, as well as the lowest RMSE and AIC values. When assessing the predictive performance of this model on the test data, an RRMSE of 23 %, a bias of −19.48 g and a MAE of 170.38 g were obtained. Furthermore, graphical analysis of the residuals for the E3 model revealed no evidence of heteroscedasticity or autocorrelation between the residuals (Fig. 6). To represent the evolution of the cumulative productions generated by the model (Fig. 6), a time (t) of 56 days from the beginning of the tapping campaign was selected. This period was chosen as a reference since it is two months after the start of the campaign and one month before the halfway point of the campaign, allowing for an early estimation of the inputs for the season. To represent the production curves, the difference between the maximum and minimum cumulative production in the time t was divided by four, and this amount was added Table 2 Base models and ADA and GADA formulations considered. Base model Parameter related to site Solution for X with initial values (t 0 , Y 0 ) Dynamic equation Bertalanffy-Richards Y=a1⋅(1−e−a2⋅t)a3 a 1 =X X0=Y0 (1−e−b2⋅t0)b3 Y=Y0⋅(1−e−b2⋅t 1−e−b2⋅t0)b3 (E1) a 2 =X X0= −ln(1− (Y0/b1)1/b3) t0 Y=b1⋅⎛ ⎝1−(1−(Y0 b1)1/b3)t/t0⎞ ⎠ b3 (E2) a 3 =X X0=ln(Y0/b1) ln(a−e−b2⋅t0)Y=b1⋅(Y0 b1) ln(1−e−b2⋅t) ln(1−e−b2⋅t0) (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 √) with L0=ln(a−e−b1⋅t0) Y=Y0⋅(1−e−b1⋅t 1−e−b1⋅t0)b2+b3/X0 (E4) Korf Y=a1⋅e−a2⋅t−a3 a 1 =X X0=Y0 e−b2⋅t−b3 0 Y=Y0⋅e−b2⋅t−b3 e−b2⋅t−b3 0 (E5) a 2 =X X0= − ln(Y0/b1)/t−b3 0 Y=a1⋅(Y0 b1)(t0 t)b3 (E6) Weibull Y= a1⋅(1−e(− a2⋅ta3)) a 2 =X X0= − ln(1− (Y0/b1) )/tb3 1 Y=b1⋅(1− (1−Y0/b1))(t/t0)b3 (E7) Hossfeld Y=a1 1+a2⋅t−a3 A 2 =X X0=t−b3 0⋅(b1/Y0−1)Y=b1 1−(1−b1 Y0)⋅(t0 t)b3 (E8) Fig. 4. Resin yields of individual grooves during the tapping season of the two methods, non-mechanized tapping method (continuous line) and mechanized circular groove method (dashed line). ´ O. L´ opez-´ Alvarez et al. Forest Ecology and Management 549 (2023) 121501 6 from the minimum to the maximum production in the period t. The represented curves are those of 187.00, 410.25, 633.50, 856.75, and 1080.00 accumulated grams at time t. Fig. 8 plots the RRMSE obtained by the model when making predictions on test data as a function of time. The model achieves an RRMSE of less than 30 % for forecasts up to 90 days (Fig. 8), which is very good for forecasts of less than one month, good for forecasts of two months and fair for forecasts of three months according to Li et al. (2013). This statistic rises above 30 % for forecasts longer than 90 days. 3.4. Mechanized circular groove method The fits to the data obtained with the mechanized circular groove method (Table 4), show favourable results but more modest compared to the non-mechanized tapping method. The R 2 values achieved with the training data range from 0.695 to 0.82. The RMSE values and AIC values range from 212.512 to 276.569 and from 16844.4 to 17496.84, respectively. When assessing the models on test data, the RRMSE was again obtained values below 30 %. The bias varied from 0.98 to 41.45 g, and the MAE ranged from 150.95 and 199.56 g. Consistently, equations derived from the Bertalanffy-Richards model demonstrate superior goodness-of-fit statistics, yielding the most optimal fit to the data. Conversely, the poorest fits are observed with the same equations as in the non-mechanized tapping method case, specifically the anamorphic and multiple asymptotic form of the BertalanffyRichards equation, as well as the Korf model by solving the “a 1 ″ parameter. The E4 or GADA form of the Bertalanffy-Richards model is found to Fig. 5. Accumulated resin yield during resin tapping season for both resin tapping methods, non-mechanized tapping method (graphic at the top) and mechanized circular groove method (graphic at the bottom). ´ O. L´ opez-´ Alvarez et al. Forest Ecology and Management 549 (2023) 121501 7 be the equation that achieves the highest level of statistical goodness of fit with the training data (Table 4). Graphical analysis of the residuals for the E4 model indicates the absence of any observable heteroscedasticity or autocorrelation among the residuals (Fig. 7). When tested against the validation data, the equation displayed the lowest RRMSE value of 29 %, a bias of 0.98 and a MAE of 150.95 g. To represent the curves of the adjustments made with the equation on the original data, it has been done in the same way as in the previous Table 3 Parameter estimates and goodness-of-fit statistics for non-mechanized tapping method. Model Train Test Parameter Estimate Standard error 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 Parameter estimates and goodness-of-fit statistics for mechanized circular groove method. Model Train Test Parameter Estimate Standard error 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 left figure was the age-dependent dynamic model E3 for non-mechanized notching and the accumulated resin yield curves of187.00, 410.25, 633.50, 856.75, and 1080.00 g at the time t of 56 days. The right figure is the Residuals vs Residuals i +1 for non-mechanized notching model E3. ´ O. L´ opez-´ Alvarez et al. Forest Ecology and Management 549 (2023) 121501 8 case. In this case, the curves represented at time t are those of 55.00, 333.75, 612.50, 891.25, and 1170.00 g (Fig. 7). The model obtained an RRMSE lower than 35 % when the prediction lag was 85 days (Fig. 8), also reducing the mean prediction error at the end of the tapping season below that obtained by the model fitted to the other extraction method. In this case, forecasts of less than one month could be described as good and those of less than two and a half months as fair (Li et al., 2013; Despotovic et al., 2016). 4. Discussion This study represents the first case of modelling individual tree accumulative resin yield throughout the resin production season. It demonstrates that methods similar to those used to model different natural growth processes can be successfully applied. The unprecedented nature of this production modelling poses a challenge to the discussion, given the lack of previous studies on resin to establish a comparative framework. Fitting models describing the trend of cumulative production during the NTFP harvest season is challenging due to the high variability in individual tree production (G´ omez-García et al., 2022; L´ opez-Alvarez et al., 2023b) and the influence that environment changes have on the yield of this renewable product (V´ azquez-Gonz´ alez et al., 2021b). Regardless of the extraction method, the equations derived from the Bertalanffy-Richards model with ADA and GADA methodologies were the ones that presented the best performance. Likewise, when comparing the statistics of the models based on the extraction methods, the model for the non-mechanized method obtained better goodness-offit than the mechanized method. The obtained results suggest that non-mechanized tapping may lead to a more pronounced yield increase in individual grooves towards the middle and end of the season when compared to the mechanized circular groove method (Fig. 4). This behavior, derived mainly from the tapping method (L´ opez-´ Alvarez et al., 2023b), accentuates the presence of an inflection point in the cumulative production curve, making it more similar to growth processes. In the non-mechanized tapping method, the equations that obtained the highest R 2 in the training phase were those derived from the Bertalanffy-Richards equation by ADA and GADA methodology (E3 and E4, respectively). The similarity in describing the trend of the data between the ADA and GADA forms of the BertalanffyRichards model is common in other forest species growth modelling works (Trim et al., 2020; Wang et al., 2020). The GADA form generally obtains better goodness-of-fit statistics, as it can generate polymorphic curves with multiple asymptotes. This statement is verified in the case of the non-mechanized extraction method, as the GADA formulation obtained the best statistics, similar to the reported by S´ anchez-Gonz´ alez Fig. 7. The left figure was the age-dependent dynamic model E4 for mechanized circular groove and the accumulated resin yield curves of 55.00, 333.75, 612.50, 891.25, and 1170.00 g at the time t of 56 days. The right figure is the Residuals vs Residuals i +1 for mechanized circular groove model E4. Fig. 8. Relative root mean square error (RRMSE) by lag of prediction for mechanized circular groove (circles) and non-mechanized notching (triangles). ´ O. L´ opez-´ Alvarez et al. Forest Ecology and Management 549 (2023) 121501 9 et al. (2008) when modeling cork growth. The goodness-of-fit statistics obtained in both fits are similar to those obtained when using the GADA form of the Bertalanffy-Richards equation to model production over time of long rotation species such as the pedunculate oak (G´ omez-García et al., 2015), and slightly lower than those obtained when modeling short rotation species (Di´ eguez-Aranda et al., 2005; Seki and Sakici, 2017) or cork production (S´ anchez-Gonz´ alez et al., 2008). The good performance of the models in terms of goodness-of-fit, as obtained in this ground-breaking study, verifies that these techniques can be used to model resin production during the season. To make these models a valid tool for modelling cumulative resin production, both in the field of activity of forest workers, in the resin tapping industry and in the resin processing industry, it would be necessary to develop this preliminary study further and to carry out a larger sample, representative of the whole environmental gradient and including all the possible casuistry under which resin pine forests can be found. Therefore, based on the above, describing the individual tree accumulated resin production of P. pinaster can be included among the uses of the von Bertalanffy model (Bertalanffy, 1949, 1957). The predictions made with the test data yielded average RRMSE values below 30 %, regardless of the extraction method. These average values can be classified as fair according to Despotovic et al. (2016) and Li et al. (2013). Simplifying and not considering the particular lag of each prediction, the models have an average accuracy of 76.81 % and 70.76 %, for the non-mechanized and mechanized methods, respectively. These accuracy percentages are slightly lower than those achieved by the growth models used in dendrometry (Prada et al., 2019) and are more similar to those obtained by S´ anchez-Gonz´ alez et al. (2008) modeling cork growth. When evaluating the accuracy percentage as a function of the lag of the predictions, it was found that a prediction made with an 80-day lag, which is usually half of the resin-tapping campaign, had an error rate of approximately 35 %. These accuracy percentages as a function of the lag are slightly lower than those achieved by dendrometric models (Gea-Izquierdo et al., 2008). The lower precision values may be due to the small sample size or the high variability of resin production over time due to dasometric and environmental variables, as data from dendrometric variables tend to show less variability in production between measurement periods. To facilitate scientific research, one potential functionality of the results obtained is their usefulness to recover missing data from the temporal weight series. Collecting intermediate weight data is a difficult task (Calama et al., 2011), and usually there are data that were not measured at the precise moment when each new strip was done but were accumulated and weighed after two or more tapping periods. In this way, the accumulated production is accounted for, but it is not known which production belongs to each period. By using this type of models, the production curves can be recreated, and the missing intermediate data can be recovered. They are particularly relevant for researchers who wish to generate and monitor production curves and have missing data in the production series. The good predictive capacity of the models makes them a valid tool for the multifunctional management of pine forests, which is feasible as long as the owners are provided with the appropriate tools to carry it out (Sabastian et al., 2019). By using these models, forest managers will be able to predict inputs before the resin harvesting season ends, allowing them to optimize the duration of the campaign to maximise the annual income from resin sales. This makes resin extraction more attractive to forest producers, creating green jobs and promoting the bioeconomy by encouraging the production of a bioproduct such as pine resin. 5. Conclusion With the increasing importance of resin production in the NTFPs and forestry sector, there is a need to further develop and expand our knowledge. This work demonstrates for the first time that age-invariant models are a valid tool for the management of pine forests under resin tapping procedures. According to R 2 , RMSE, AIC, RRMSE, Bias and MAE, the ADA and GADA formulations of the Bertalanffy-Richards model performed best in modelling cumulative resin production. This implies the inclusion of accumulative resin yield of individual trees as one of the many phenomena that can be modelled with these equations. Furthermore, the ability to accurately predict final production once the campaign has started using these models represents a significant advance in multifunctional forest management. In addition, this type of model applied to cumulative annual resin production can provide a basis for reconstructing missing data series with a high degree of accuracy. To further advance this sector, public or private initiatives should provide support and follow the example of forest growth models by establishing a network of permanent plots for parametrising national-scale resin production variables and developing more robust models. Enabling advancement in the effective joint management of wood resources and resin production, facilitating the creation of technical–economic models and the implementation of innovative, cutting-edge methodologies in constructing said models. Founding and acknowledgements This work was supported by the Spanish Government (“ACREMA”, MAPA/AEI-Agri/FEADER, UE) [O00000226e2000043659], the Galician Government (Xunta de Galicia) with a grant for Competitive Reference Groups ED431C-2021-27 and the pre-doctoral contract Campus Terra-USC 2023. The authors thank FORESIN and CIF Lourizan, for the field assessments. CRediT authorship contribution statement ´ Oscar L´ opez-´ Alvarez: Conceptualization, Methodology, Formal analysis, Writing – original draft, Visualization. Luis Franco-V´ azquez: Methodology, Writing – review & editing. Manuel Marey-Perez: Conceptualization, Methodology, Formal analysis, Resources, Writing – review & editing, Supervision, Project administration, Funding acquisition. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability Permission to share the data must be obtained from the institutions that funded the research References Bailey, R.L., Clutter, J.L., 1974. Base-age invariant polymorphic site curves. For. Sci. 20, 155–159. https://doi.org/10.1093/forestscience/20.2.155. Bertalanffy, L.V., 1957. Quantitative laws in metabolism and growth. Q. Rev. Biol. 32, 217–231. https://doi.org/10.1086/401873. Caglar, M.U., Teufel, A.I., Wilke, C.O., 2018. Sicegar: R package for sigmoidal and double-sigmoidal curve fitting. PeerJ 6, e4251. Calama R., Miina, J., de-Miguel, S., Bonet, J. A., Mounir, F., Tom´ e, M., MartínezJaúregui, M., Herruzo, C., Peltola, R., Salo, K., Kurttila M., Hern´ andez-Rodríguez, M., Martín-Pinto, P., S´ anchez-Gonz´ alez, M., 2020. Data & models: importance of assessing and forecasting non-wood forest products in Europe, in: Vacik, H., Hale, M., Spiecker, H., Pettenella, D. Tom´ e, M., Non-Wood Forest Products in Europe. Ecology and management of mushrooms, tree products, understory plants and animal products. Outcomes of the COST Action FP1203 on European NWFPs, BoD, Norderstedt, pp. 43-78. Calama Sainz, R., Tome, M., S´ anchez-Gonz´ alez, M., Miina., J., Spanos, K., Palahi, M., 2011. Modelling Non-Wood Forest Products in Europe: a review. Forest Syst. 3 (4), 69. ´ O. L´ opez-´ Alvarez et al.