A tree- and climate-dependent growth model to predict mature annual cork thickness under different climate change scenarios

Climatic factors drive the annual growth of cork and the subsequent increase in its thickness, which, in addition to porosity, determines the price of cork. Therefore, the simulation of cork thickness is a crucial module of forest growth simulators for cork oak stands. As the existing cork growth models are independent of climatic factors, cork thickness under different climate change scenarios could not be simulated using these models. The primary objective of this study was to develop a climate-dependent tree model to predict annual cork growth. We also verified the hypothesis that the effects of climate change on cork annual growth are nonlinear, and vary with the cork age and thickness. Due to the limited amount of work developed around this topic, we evaluated three candidate models and selected the one that presented best prediction performance as the base model. A set of climate variables that characterized annual climatic conditions were tested in the base model parameters. The resulting climate-dependent model was referred to as the fixed-effects model, and used to initialize a mixed-effect model which accounted for the nested structure of the data. We considered two random effects—the plot and the trees inside the plot. Annual precipitation and the Lang index (ratio between annual precipitation and mean annual temperature) were the variables that showed best results when included in the model parameters. Using a ratio of the variable to cork thickness recorded during the previous year, in both cases, suggested a decline of the positive effect of annual precipitation and the Lang index for increasing cork thickness. The models developed in this study predicted the cork thickness of individual trees based on the cork age and under different climate change scenarios. Therefore, they can be used in forest growth simulators for forest management and research purposes.


Introduction
According to the FAO, non-wood forest products are 'products of biological origin other than wood, derived from forests, other wooded land and trees outside forests (FAO 1999). The Mediterranean agrosilvopastoral ecosystems-Montado and the Dehesa-provide several non-wood forest products, with the bark of Quercus suber L. (the cork) being one of the most important products originating from the tree component of the ecosystem (Moreno et al. 2017;Vacik et al. 2020). The predicted increase in the severity and frequency of drought events in the Mediterranean region owing to climate change is a major concern for forest managers and owners, farmers, and ecologists (e.g., Soares et al. 2015). The impact of these events on the Montado and Dehesa ecosystems, in terms of species biodiversity, ecosystem resilience, and socioeconomic indicators, has been well documented (e.g., Costa et al. 2011;Jongen et al. 2013;Iglesias et al. 2016;Oliveira et al. 2016).
The effects of climatic factors on cork growth, namely precipitation and temperature, were first reported by Caritat et al. (1996Caritat et al. ( , 2000. Further and more recent studies suggested the role of a wide range of edaphoclimatic variables on cork growth, such as seasonal climate conditions (e.g., Costa et al. 2002;Sánchez-González et al. 2007;Pizzurro et al. 2010;Oliveira et al. 2016;Paulo et al. , 2021Ghalem et al. 2018), groundwater availability (Mendes et al. 2016) or understory management (Faias et al. 2018(Faias et al. , 2019. Therefore, climatic factors drive annual cork growth, reducing average cork growth in the Mediterranean region, particularly under severe drought conditions (Oliveira et al. 2016). Although climate change does not alter the chemical constitution of cork (Leite et al. 2020), it does affects the cork rings density (Costa et al. 2022) and the annual growth rate (e.g., Caritat et al. 1996Caritat et al. , 2000. These, in turn, affects the caliber (or thickness) of the extracted material, with impacts to the value chain by a reduction of the percentage of usable material for the production of cork stoppers (Paulo et al. 2021), which ultimately influences cork price and economic return for the farmers .
Despite being harvested for centuries, non-wood forest products are seldom integrated into forest management plans. Therefore, investigations and decision support tools development are necessary to ensure the long-term sustainability of these resources and realize their economic potential (Chamberlain et al. 2019). In this regard, forest models and simulators, particularly those that are consistently updated to incorporate recent scientific findings regarding forest ecosystems response to external drivers (climate and hazards), are crucial (Peng et al. 2006).
There have been several studies on cork since the late nineteenth century, with results and observations first summarized in the middle of the twentieth century by Natividade (1950). Over the last four decades, numerous studies have focused on the composition and growth of cork, which has been summarized by Pereira (2007). Nonetheless, only a few forest simulators have been available as support tools for the management of Montado and other areas with the prevalence of cork oak.
Till date, only three models have been developed: Suber v5.0 (Paulo 2011), Corkfits (Ribeiro and Surový 2011), and Alcornoque (Sánchez-González et al. 2015). Each of these models simulates annual growth and thickness using various growth functions, fitted with distinct datasets and statistical approaches, such as Almeida et al. (2010) for Suber, Ribeiro et al. (2003) for Corkfits, and Sánchez-González et al. (2008) for Alcornoque. These forest simulators, or a few of the forest models that they include, have been increasingly used for studying the management of the Montado ecosystem (Borges et al. 1997), optimization of cork debarking schedule (Palma et al. 2015;Pasalodos-Tato et al. 2018), and assessment of ecosystem services Rodrigues et al. 2021). Despite the differences between the existing models, all of them estimate cork growth and thickness independent of climatic conditions. Therefore, these models cannot provide accurate estimates under different climatic scenarios. To the best of the knowledge of the authors, no cork growth model of cork oak forest simulator incorporates climate variables for the prediction of cork growth. Therefore, the objective of the present study was to develop the first treelevel annual cork growth model that included model parameters depending on climate variables. The addressal of this objective also allowed to answer the following questions: (i) What are the optimum annual climatic conditions for annual cork growth, (ii) Are these conditions dependent on the cork age and thickness?

Materials
Cork samples were collected from trees in 35 permanent monitoring plots that were installed in even-or unevenaged rain-fed cork oak stands distributed throughout the Portuguese cork production area. The stands and thus the plots were managed using traditional forest or silvopastoral management practices, suggesting that the stands were never irrigated. Furthermore, understory composition and management varied among stands, with understory control being made by cattle grazing (cow or sheep) or mechanical equipment mechanically (bush cutter or disc harrow) at a frequency of 4-15 years. Details of these plots can be obtained from Paulo et al. (2021).
Cork samples (20 × 20 cm 2 ) were collected from each tree at a height of 1.3 m. During the first (1996)(1997)(1998) and second (2005-2010) sampling periods, a total of 645 and 642 cork samples were collected from the stands, respectively.
After sampling, the cork samples were processed in the laboratory according to standard post-harvest procedures. Briefly, the samples were immersed in boiling water for 1 h at 100 °C and atmospheric pressure, followed by drying at an ambient temperature. Details of the cork post-harvest process can be found in Pereira (2007). Thereafter, transverse sections of each sample were polished using a belt sander (FAI™ Model LCU) fitted with a p180 sandpaper band (Norton™). This facilitated the posterior identification of annual rings, carried out the following procedure described by Paulo et al. (2021). The two half cork growth rings, corresponding to the years of debarking and located at the two extremes of the samples, were not considered in this study. This implied that cork samples collected at an age of t years comprised t − 1 complete cork growth rings.
The final dataset included 12,351 measurements of complete annual cork rings width or increment (ict), out of which only 11,208 were used in this study. The remaining 1143 annual growth measurements (9%) were discarded because of the difficulties in precisely identifying the limits of cork rings. Each cork ring was associated with the cork age at the growth year (t). The distribution of the number of measurements per cork ring age was similar for rings of ages 1 to 8 years, i.e., 11% observations in the dataset for each age.
However, a reduction in the number of measurements was observed for rings of ages 9 or more years. This could be attributed to the fact that the majority of the stands were debarked when the cork exhibited 9 years of growth, thereby including only eight complete annual rings. The dataset included 6%, 3%, and 3% measurements corresponding to the 9th, 10th, and 11th-19th complete growth rings, respectively, and the remaining growth rings corresponded to ages 18 and 19 years. Similar to previous studies, ict range of values decreased with increasing t (Fig. 1).
For each year and plot included in the dataset, we obtained data on annual precipitation (Pt t ) and mean annual temperature (Tm t ) from the year t to determine annual climatic conditions (Table 1). Annual variables were computed from October of the year t − 1 to September of the year t. The Lang index (LI), computed as the ratio between annual precipitation and mean annual temperature (Lang 1915), was used to classify climatic conditions into arid (LI < 40), semi-arid (40 < LI < 60, humid (60 < LI < 100), and perhumid (LI > 100). The dataset, thus, included varying annual climatic conditions, with values of LI ranging from 9.8 (arid) to 68.2 (humid) and a mean value of 42.3. Accordingly, 49.5%, 43.4%, and 7.1% of the years were classified as arid, semi-arid, and humid, respectively.
Two additional hypotheses were formulated and verified, which resulted in the identification of new variables as functions of previous ones. The first hypothesis tested if the impact of annual climatic conditions on cork growth was linear or presented an optimum value. To verify this hypothesis, climate variables were expressed in linear and quadratic forms (e.g., a 0 Pt 2 t + a 1 Pt t + a 2 ). The second hypothesis stated that the effects of annual climatic conditions on annual cork growth were dependent on cork age (t) and cork thickness (c t ). For this purpose, the ratios between linear or quadratic annual climatic variables and cork age ( Pt t t or a 0 Pt 2 t +a 1 Pt t +a 2 t ) or cork thickness at the year t ( Pt t ct t or a 0 Pt 2 t +a 1 Pt t +a 2 c t ) were considered.

Selection of the base model
The model, namely the individual tree annual cork growth model, was developed as a differential equation model that estimated the accumulated cork thickness at the age of t + 1 years (ct t+1 ) as a function of the accumulated cork thickness at age of t years (ct t ) and cork age. The model parameters were expressed as functions of climatic variables. For details of the methodology employed for the formulation of forest growth models using differential equations, see Burkhart and Tomé (2012). The candidate models were selected based on the findings of Almeida et al. (2010), who developed a system of equations to predict the evolution of mature cork caliber of an individual tree over time using measurements of the first cork samples collected between 1996 and 1998. After comparing 13 models [see Table 4 of Almeida et al. (2010)], the authors selected the algebraic differential equation (ADA) log-logistic model with parameter 'a' representing the tree-specific parameter (LL_a; Eq. 1). Similar results were obtained using two other ADA models with one treespecific parameter: the log-logistic model with parameter 'b' representing the tree-specific parameter (LL_b: Eq. 2) and Lundqvist model (L_a; Eq. 3) (Lundqvist 1957) with 'a' representing the tree-specific parameter. Considering the lack of studies on this theme, these three models were used as candidate models in the present study: Each of the three candidate models was initially fitted with all parameters independent of climatic variables using the nonlinear ordinary least-squares method in the MODEL procedure of the SAS software (SAS Institute Inc. 2014). This allowed the establishment of initialization values of parameters that were only available for the LL_a model (Almeida et al. 2010). The best-fit model, considering mean square error and adjusted r 2 , was selected as the base model.

Fitting of the fixed-effects model
The next step focused on the refitting of the base model now including one or two model parameters expressed as a function of climatic variables that characterized climate growth conditions in the year t + 1. This was performed using the nonlinear ordinary least-squares method implemented in the MODEL procedure of the SAS software (SAS Institute Inc. 2014).
Models that exhibit high collinearity must be avoided as they can lead to unstable estimations due to over-parameterisation. For this reason, each climatic variable was separately introduced into one of the model parameters, and a collinearity diagnostics was carried out. This resulted in the fitting of 48 models. The models with a condition number above 30 were considered to exhibit collinearity and discarded (Myers 1990). For models with a condition number equal or under 30, the significance of all coefficients was assessed by asymptotic t tests. Variables that were significant upon addition to one parameter were then combined with the two parameters of the base model. To access the improvement from adding climate variables to the base model parameters, a relative change metric was calculated for the mean square error, formulated as where MSE_0 is the mean square error of the base model, and MSE_c is the mean square error of the each fitted model after expressing one or two model parameters as a function of climate. The smaller the %MSE_0 value, the better the candidate model performance.
The fitted models were then compared using fitting and prediction statistics, considering the press residual (r p ) (Myers 1990). Fitting statistics included mean square error (MSE) and adjusted r 2 (adj-r 2 ), whereas prediction performance was evaluated using mean press residual (Mrp) for bias assessment and mean of the absolute value of the press residual (Marp) for precision assessment. The proportion of variation explained by the model, which is designated by model efficiency (ef), was estimated as follows: where r pi denotes the press residual for the observation i, y i denotes the value of the observation i, y denotes the average value of the observations, and n represents the number of observations.
Owing to the importance of accessing model performance across different cork ages and range of annual climatic conditions, the prediction statistics were also plotted across the range of cork age and climatic variable values represented in the data set, with a focus on annual precipitation and LI classes. Thereafter, the homoscedasticity and normality hypotheses were assessed by examining the plot of studentized residuals against the predicted values and by analyzing the qq plot, respectively. The model with the best fitting and prediction behavior was used as the final fixedeffects model.

Fitting of the mixed-effects model
The data used in this study were characterized by a nested structure, because of the presence of several trees inside plots. Compared to Sánchez-González et al. (2021), the present study did not consider the region random effect, as the definition of cork regions did not apply to Portugal (Pereira 1999). The mixed model methodology (Pinheiro and Bates 2000;Fang and Bailey 2001) has been successfully employed for modeling cork growth and thickness using data sets of these characteristics (Sánchez-González et al. 2007. Therefore, the final fixed-effects model, developed in "Fitting of the fixed-effects model" section, was refitted after it was redefined as a mixed model. For this purpose, the model parameters were reformulated as mixed parameters by including one or two random effects in each model parameter, separately and simultaneously. The two random parameters represented the effects of plots and the trees inside the plots. A total of nine mixed models were fitted and compared using the SAS macro NLINMIX (SAS Institute Inc. 2004). The option of expanding on the random effects at their current empirical best linear unbiased predictors (EBLUP) was used. For a few models, convergence implied an increase in the number of iterations, which was performed by including option 'MAXIT = 50'. Basic assumptions for the nonlinear mixed model fitting included the multivariate normal distribution of random effects and the error term vector. When the plot of the resulting studentized residuals versus the predicted values suggested heterogeneous intra-individual error variance, variance functions were verified using the methodology described in Paulo et al. (2011).
The resulting mixed models were compared with each other and with the final fixed-effects model. Due to their nested structure, resulting models were compared by likelihood ratio tests (LRT), Akaike's information criterion value (AIC), and the percentage reduction in AIC compared to the fixed-effects model (%AIC_0). Considering the latest, the smaller the %AIC_0 value the better.
We developed two alternative models: one was restricted to fixed-effects parameters and the second was a mixed-effects model. The former may be used when the calibration of the mixed-effects model is not possible as, under these circumstances, predictions obtained from models fitted using OLS are less biased than the ones resulting from mixed model estimations under population-specific responses. This may occur as in addition to collecting cork samples during forest inventory, investing resources for the preparation, boiling and cork ring measurement are required for calibrating the mixed-effect model. Therefore, the two alternative models account for different use of the results of the present study.

Annual cork thickness prediction under different climate change scenarios
To demonstrate the simulation of the impact of different climate scenarios on cork thickness, the fixed-effects model was used to simulate cork caliber increase over distinct 13-year simulation periods. Cork growth was estimated considering the values of annual precipitation and average temperatures observed at the period 1971-2003, and the simulated values obtained from the RCP4.5 scenario for the period 2011-2043. Climate data were obtained from the Portal do clima (http:// porta ldocl ima. pt/ en/).
A second simulation was carried out that aimed to demonstrate the impact of a year with severe droughts occurring on different cork ages. Five series of annual cork growth increment were plotted using the simulated climate data (RCP4.5 scenario) for the period 2011-2023 (cork age 1 to cork age 13), modified to include the effects of a severe drought year on the 2nd, 3rd, 5th, 7th, and 9th years of cork age (Fig. 2). The values used for the characterization of the severe drought year were recorded in 2005 in Coruche, Portugal: annual precipitation of 329 mm and an average temperature of 15.1 °C. For these simulations, an initial cork thickness value of 3.9 mm at 1 year of age was considered. Fig. 2 Annual precipitation (mm) and mean annual temperature (°C) provided by the RCP4.5 climate scenario (http:// porta ldocl ima. pt/ en/) for the period 2011-2023, corresponding in the model to cork age 1 to cork age 13. *Identifies the years whose values were modified, each one separatly for each simulation, to simulate a severe drought occurrence. The modifications implied considering an annual precipitation of 329 mm and an average temperature of 15.1 °C, instead of the original value from the RCP4.5 climate scenario

Selection of the base model
The statistics resulting from the fitting of the three candidate models suggested the exclusion of the LL_b model. Considerable differences were not reported for the other two models-LL_a and L_a (Table 2).
The LL_a model was selected as the base model (Table 3) as it presented the 5th and 95th press residual percentile values close to zero. Parameter estimates of the three candidate models are presented in "Appendix".

Fitting of the fixed-effects model
The fitting of the LL_a model with one parameter as a function of a climatic variable (48 models fitted) resulted in 15 alternative models (Table 4), 8 with the 'b' parameter and 7 with the 'c' parameter dependent on annual climatic conditions. The remaining were discarded because of collinearity or the non-significance of one or more parameters. Moreover, only one of the alternative models included the mean annual temperature, whereas the remaining 14 models exhibited parameters as a function of annual precipitation or LI. We also observed that for most of the alternative models, precipitation and LI were expressed in parabolic forms and/ or as their ratios to cork thickness or age. The combination of climatic variables, previously included in one model parameter, simultaneously in two parameters resulted in the fitting of 71 models. Among these, 29 resulted in alternative models with all significant parameters and no collinearity ( Table 4).
The observation of the models from Table 4 showed that the variable mean annual temperature exhibited minimum contribution compared to other variables, as it was included in only 7 out of 44 models, and only in the 'b' growth parameter. On the other hand, only one model did not include the annual precipitation or the Lang index. The importance of the 'c' parameter (inflection point) for the cork annual growth modeling was also demonstrated in the first 15 models of Table 4, as all included this parameter expressed as a function of climate.
Comparing candidate models based on fitting and validation statistics and after sorted for increasing values of %MSE_0, revealed that the model with the 'b' parameter as a function of the ratio of LI and cork thickness and the 'c' parameter as a function of the ratio of annual precipitation and cork thickness was selected as the fixed-effects model ( Table 5). In fact, when compared with the base model, the model showed a reduction of MSE (%MSE_0 = 80.8%).
It was also noticed that the models ranked from second to sixth in Table 4, included the 'c' parameter which expressed an upward parabolic function of the ratio between the LI and the ct. Parameter estimates of the second and third best ranked models are presented in "Appendix".
The plot of the studentized residuals versus the predicted values did not evidenced within-individual error heterogeneous variance. Furthermore, no differences were observed in the plots of studentized residuals versus cork age, LI or annual precipitation (Fig. 3). The normality of the residual hypothesis was also verified using the qq plot.

Fitting of the mixed-effects model
The LL_a fixed-effects model was refitted after reformulation as a mixed-effects model. In this step, nine mixedeffects models, with a different number of mixed and random parameters, were fitted and compared. We did not observe a considerable decrease in AIC value for all models, particularly the models that included the 'b' parameter as a single random parameter (Table 6). However, reduced AIC values were observed when both the 'b' and 'c' parameters were defined as mixed parameters. Table 2 Model fitting and validation statistics of the three candidate models MSE: mean square error; adj-r 2 : adjusted r 2 ; P 5 : 5th percentile of press residuals; Mrp: mean of the press residuals; P 95 : 95th percentile of press residuals; Marp: mean of the absolute press residuals; ef: model efficiency; LL_a: log-logistic model with 'a' defined as tree specific; LL_b: log-logistic model with 'b' defined as tree specific; L_a: Lundqvist model with 'a' defined as tree specific. * selected as base model  A reduced AIC value was associated with the selected mixed-effects model selected (rank = 1 in Table 6), which exhibited 'b' and 'c' parameters as mixed parameters while incorporating both plots and trees inside the plots as random effects (Table 7).  Table 2); MSE: mean square error; adj-r 2 : adjusted r 2 ; P 5 : 5th percentile of press residuals; Mrp: mean of the press residuals; P 95 : 95th percentile of press residuals; Marp: mean of the absolute press residuals; ef: model efficiency. Models with non-significant parameter estimates or a condition number > 30 are not presented in the table *Indicates the selected model The plot of the studentized residuals versus the predicted values for this model did not reveal heterogeneous intra-individual error variance. We observed only slight differences between this plot and that of the fixed-effects model. Additionally, no trend was observed while plotting the studentized residuals versus the observed values of cork age, annual precipitation, mean annual temperature, and LI (Fig. 4). The hypothesis of the normality of the residuals was also verified using the qq plot. Figure 5 exhibited three plots, considering the simulation of cork annual growth of three distinct trees. The trees are characterized by distinct cork thickness values at 1 year of age. Each plot includes the simulation of cork thickness increase for the following 13-year periods: 1971-1983, 1981-1993, 1991-2003, 2011-2023, 2021-2033, and 2031-2043 (each 13-year period is represented by one series).

Annual cork thickness prediction under different climate change scenarios
Climate change impact was observed in two aspects. From one hand, the increase of cork age that allows reaching a value of 30.7 mm (horizontal line in Fig. 5), considered as the minimum value, after boiling, for the production of natural cork stoppers. For trees characterized by a smaller value of initial cork growth, the value was already 11 years for simulations carried out with 1970-1983 climate. In future climate scenarios, this value increased even more for 13 years (Fig. 5 left). For trees characterized by median cork thickness values at 1 year of age, an increase from 9 to 10 years is already simulated when considering moving from the 1970-1983 period to the 1991-2003 period. A delay of 11 years was simulated when considering the RCP4.5 climate scenarios (Fig. 5 middle). For trees characterized by larger values of the model initialization, values under 9 years were simulated, the minimum valued allowed by the Portuguese legislation (Fig. 5 right).
Climate change impact was also observed by the considerable reduction in cork thickness estimated for the same cork age, when these were carried out using distinct climate scenarios. Upon estimating the percentage of cork thickness reduction, at 9 years of age, and considering the 1971-1983 period as a base value, the reduction reached values ranging from − 9.9 to − 13.2% under climate that has already characterized the period 2011-2023. Under RCP4.5 climate scenarios for 2031-2043, estimated values show that a reduction of cork thickness, at 9 years of age, is expected to reach − 10.5 to − 15.2%. In all simulations, trees characterized by lower values of model initialization presented the larger values of percentage of cork thickness reduction. Figure 6 demonstrated the impact of a year characterized by a severe drought in annual cork growth rates, when this occurred for different cork ages. When it was simulated for cork age of 2 or 3 years, the estimated annual cork growth showed a clear reduction, one that is more smoothed when the severe drought year occurs latter in the cork growth cycle (cork age 7 or 9 in Fig. 6). These values, outcome of the inclusion in the model parameters of climate variables, expressed as the ratio of climate variables and cork thickness.

Discussion
This study draws upon cork growth models developed by Almeida and Tomé (2008), Almeida et al. (2010), Sánchez-González et al. (2007;, to develop the first individual tree-and climate-dependent annual cork growth model. The range of annual precipitation and temperature included in the present data set will allow the developed model to be applied across the natural distribution areas of cork oak species. Owing to its characteristics, annual cork growth was modeled using differential equations (Burkhart and Tomé 2012), where a variable measurement taken at one time is used to Simulation of cork thickness increase with the increasing cork age obtained using the proposed fixed-effects model. From left to right, the three plots show the simulation considering a tree with cork thickness values at 1 year of age equal to 2.8, 3.9, and 5.2 mm, corresponding to the 1st, 2nd, and 3rd quartile values obtained from the data set respectively. On each plot, six series were included, obtained considering the observed values of annual precipitation and average temperature for the periods 1971-1980, 1981-1990, and 1991-2000, and the simulated values (the RCP4.5 scenario) for the periods 2011-2020, 2021-2030, and 2031-2040. Horizontal line indicates the minimum value of cork thickness, after boiling, suitable for natural cork stopper production (30.7 mm) predict future or past variable estimates. Based on the documented relationship between annual cork growth and annual climatic conditions, the primary objective of the present study was to model cork thickness accumulated in the year t + 1 based on the thickness in the previous year t and the climatic conditions of the growth year t + 1. The annual time step adopted for model formulation would be helpful to forest growth simulators and forest managers, whereas producers and users of simulated data would find the model suitable for sensitivity analysis or optimisation methodologies focused on the establishment of the debarking rotation period under climate change scenarios, an objective ahead of that carried out by Palma et al. (2015) and .
The present model was based on the comparison of three alternative algebraic difference equations (ADA). We used this methodology instead of using dynamic models derived from the generalized algebraic differential approach (GADA) (Cieszewski 1994;Cieszewski and Bailey 2000) because of the following reasons: (1) Almeida et al. (2010) obtained similar results using both ADA-and GADA-based models; (2) Sánchez- González et al. (2007) reported that five of the candidate GADA models could not converge in parameter estimates, even with an increased number of maximum iterations, reduced convergence criteria, and altered initial parameter values, suggesting that 'both the dynamic models derived using ADA and GADA adequately described the data'; and (3) the parsimony principle and its contribution to future applications of the model in forest management and forest simulators such as the SUBER v5.0 model (Paulo 2011).
There is an ongoing debate on the effects of increased water scarcity under climate change scenarios, particularly affecting the Mediterranean region, on the farming and irrigation practices employed (Iglesias et al. 2007(Iglesias et al. , 2011Tramblay et al. 2020). However, few studies have investigated the impact of irrigation on the annual growth in cork. Despite the values reported by Poeiras (2022) and Poeiras et al. (2022) for the sum of the annual precipitation and annual irrigation do not exceed 738 mm, the fact that the irrigation was carried out during springs and summer periods must not be overlooked. This irrigation pattern has resulted in a shift in the pattern of seasonal water availability characteristic from Mediterranean climates, and for of the data set used in the present work. Support for this concern was provided by Paulo et al. (2021) that revealed the effects of seasonal precipitation, and Oliveira et al. (2016) that highlighted the high sensitivity of annual cork growth to later spring precipitation that Poeiras (2022) has amplified with irrigation. As a result, the application of this model should be restricted to trees in rain-fed land, as the data set included only these plot characteristics.
The model developed in this study captured, for the first time in a mathematical growth model, the climate effect on cork growth (e.g., Caritat et al. (1996Caritat et al. ( , 2000, Costa et al. 2002, Sánchez-González et al. 2007, Pizzurro et al. 2010, Oliveira et al. 2016, Ghalem et al. 2018, Costa et al. 2022. It confirmed, while applied for cork growth simulation, that when comparing cork growth under climate conditions from the last decades of the XX century, with cork growth under current climate conditions or future climate scenarios under climate change, a reduction of cork thickness is observed. It also showed that the maximum impact of drought years occurs when total cork thickness is thinner. Despite the importance of cork age in the response to a severe drought year was incorporated in several alternative models, the accumulated cork thickness at a given age showed the best prediction results and exceeded the first. This shows, for the first time, that accumulated cork thickness could be a more suitable variable for the determination of the suitable conditions for debarking, in comparison to cork age that is currently used by the national regulation existing in cork producer countries. This work confirmed that under climate change scenarios, if aiming to maximize the amount of cork suitable for the production of natural cork stoppers, the debarking rotation period must be extended (Leite et al. 2019). Alternatively, landowners and forest managers might assent the harvest of an increasing amount of cork suitable for the production of other cork products such as cork discs, even if these are usages associated to lower price categories. The model developed allows the simulation of the cork thickness evolution for distinct climate scenarios. Under the availability of cork forest inventory data, characterizing a homogeneous management unit, the usage of the proposed model allows to account and optimize forest management decisions related Fig. 6 Annual cork growth increment with increasing cork age. Each series includes one severe dry year at different cork age. The initial cork thickness at 1 year of age was 3.9 mm (median value obtained from the data set) to cork debarking calendars, a step forward for an area of research still seldom addressed (Borges et al. 1997;Palma et al. 2015;Pasalodos-Tato et al. 2018).

Conclusion
For the first time, climate effect was included in a mathematical annual cork growth model. This was achieved by the development of an algebraic differential model, based on the log-logistic model with parameter 'a' expressed as a tree-specific parameter, and with both the growth ('b') and inflection ('c') parameters expressed as a function of annual climate variables.
The models proposed in this study included the impact of annual precipitation and mean annual temperature, by their ratio expressed using LI. They were also able to demonstrate that climate effects, positive or negative, are less impacting for increasing values of total cork thickness.
The usage of the model for the simulation of cork annual growth allowed to demonstrate and quantify the effect of distinct climate scenarios, including simulated climate change scenarios. The occurrence of severe drought, during the first years after debarking, was demonstrated to have a stronger effect in the reduction of annual cork growth. Under climate change, and aiming at maintaining the production of cork suitable for the production of natural cork stoppers, farmers may have to delay the debarking and extend the debarking rotation period.
Our model results in an important tool for forest managers, and may be included in forest simulators.

Conflict of interest The authors have no conflict of interest to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.