Empirical Modelling of Stem Cambium Heating Caused by Prescribed Burning in Mediterranean Pine Forest

: Little is known about the interactions between the variables involved in the post-ﬁre response of Mediterranean pine species to prescribed burning (PB). Thus, it is essential to develop an empirical model in order to assess the inﬂuence of tree and stand attributes, burn season, and ﬁre severity on the probability of stem cambium damage occurring. Prescribed burnings were conducted in different seasons and areas covering a wide climatic and ecological range. Potential explanatory variables were measured. A random effects hurdle model framework was used to evaluate the temperature duration above 60 ◦ C as a proxy for stem cambium damage at tree scale. The results showed signiﬁcant differences in cambium damage between the PB seasons. Pinus nigra was more resistant than other pine species. Bark thickness was critical for protecting cambium. Volume of crown scorch, percentage of stem scorch, and maximum outer bark temperature were directly related to temperature duration above 60 ◦ C in the cambium. Prescribed burning conducted under tree canopy in Mediterranean pine species generally results in a low level of cambium damage. Empirical models could help managers to predict the effects of PB and thus select the most suitable prescriptions.


Introduction
The number of catastrophic forest fires is growing at an alarming rate worldwide [1].About 85% of the total burned area in Europe is located in the Mediterranean region [2], which is particularly prone to wildfire.In this context, prescribed burning (PB) is used as a management tool, in combination with other landscape-scale vegetation interventions, to modify fuel load and connectivity, and limit the extent of high-severity wildfires [3,4].Prescribed burning is also used to emulate the regime of light surface fires, historically frequent in some of these ecosystems [5,6].
Despite its benefits, PB can damage the crowns, roots and stems of trees [7][8][9].Stem damage occurs when phloem and cambium tissues are lethally affected by heat [10][11][12][13].This can lead directly to tree death, especially in seedlings and saplings [14,15], and can also contribute to delayed mortality, often in combination with crown or root damage [16,17].Stem damage is central to the two leading hypotheses used to explain the immediate impact of fire on trees and also second-order effects and mortality: cambium necrosis and xylem hydraulic dysfunction [18][19][20].Field experiments have shown that non-lethal stem injury can affect sap flux [21], stomatal conductance and water use efficiency [22,23].Temporary alterations in secondary metabolites [24] involved in tissue restoration or defence against insect attack [25] have also been observed in association with stem damage.
The physical processes involved in the formation and uneven location of fire scars on the tree bole (e.g., [26,27]), particularly regarding the role played by the bark (e.g., [12,[28][29][30][31][32]), have received wide attention.In addition, a number of models have been developed in order to estimate heat transfer to the cambium and phloem and the possible necrosis of these tissues (e.g., [33][34][35]).Empirical modelling of the probability of fire-induced stem injury is less common than physical modelling [8,[36][37][38][39][40] and has mainly focused on hardwoods [41].The empirical approach can complement other methods and can be useful for decision making regarding PB.Forest managers are often faced with the challenge of determining suitable burning prescriptions to ensure that the multiple objectives of PB are achieved while minimizing the negative environmental impacts and maximizing the resilience of the treated stands.Furthermore, predicting the level of stem scarring caused by PB is of economic significance, due to the lower value of damaged timber [42][43][44].Bio-physical approaches have increased the understanding of fire-induced stem damage and the associated tree response, but they have also revealed large gaps in the knowledge (e.g., [17,20,[45][46][47][48][49][50][51]).
In the above context, this study explores an empirical approach to determine which variables affect cambium damage in tree stems during PB, using the duration of temperature above 60 • C (tC 60 ) in the cambium as the response variable.Although exposure to this temperature is generally considered to lead to cambial necrosis [34,52,53], exposure to lower temperatures can also be lethal to cambium [54].Furthermore, occasional heating of the inner bark alone may not induce a response in the cambium if it is not prolonged [7,55], and thus may be a poor indicator of damage if the duration of the heating is not considered [56,57].
As the model is intended for use in PB planning, it is potentially useful to include variables that can be controlled by forest managers and that can impact both fire behaviour and effects.Specifically, burn season has a significant influence on the tree physiological condition and its response to fire [58][59][60], although the effect on post-fire tree mortality prediction seems inconsistent [16,61,62].Seasonal processes such as fine root production and carbohydrate storage [63], foliage vulnerability [60], and water shortage [64] can act independently or can interact following burning.Moreover, some tree variables may be particularly useful for predicting the likelihood of stem damage.The overwhelming importance of bark thickness in protecting cambial tissues against fire (see references above) leads to this variable being the first choice for predicting the probability of stem damage, although the closely related stem diameter is generally more useful for determining the likelihood of stem scarring [37,65].Stand attributes, such as species, ecotype composition and stand density, can strongly influence post-fire survival [66][67][68][69][70][71] and represent potentially important candidates to be explored.Stand density affects tree competition, which acts on pre-fire vigour, frequently affecting post-fire responses and survival [72][73][74].In addition, stand structure and the associated fuels are related to potential fire intensity [75], fire severity [76,77] and post-fire mortality [78].
Some widely used surrogates of fire damage in the stem and crown may also be useful for predicting the likelihood of damage.Stem char height and percentage of bole char height have commonly been used as a proxies for potential cambium damage in post-fire mortality prediction [57, [79][80][81][82].However, the use of these variables has been criticized due to their frequent poor or inconsistent performance [17,83,84].Furthermore, stem char height has also been used as a proxy for flame length and fireline intensity [85].Bark char severity [86], an indicator of char depth used in mortality prediction studies [16,57,70], has been found to be significantly related to the probability of cambium necrosis [8].Crown scorch height and percentage of crown volume or length scorched are by far the most commonly used variables to predict post-fire tree mortality [83,87,88], both alone and in combination with stem char height [62,89].Since crown scorch does not describe the actual amount of crown foliage or buds affected by fire [90], monitoring litterfall biomass [31,91] can also help to determine the level of crown damage and can serve as a proxy for fire severity.
Accordingly, we conducted a series of prescribed fire experiments in different Mediterranean pine forest ecosystems and hypothesized that tree cambium heating and the damage caused by a spreading fire may be affected by three types of variables: tree and stand variables; fire severity and damage proxies; and prescription-linked variables such as burning season.The specific research objectives are as follows: (i) to ascertain the influence of tree and stand variables that protect stem cambium from fire-caused damage; (ii) to explore the potential capacity of fire severity proxies to reflect cambium heating during PB; and (iii) to assess the effect of burn season on cambium heating.For this purpose, we propose the use, for the first time in this context, of a random effects hurdle model framework.

Study Sites
Data were collected from 5 locations in the Iberian Peninsula (Figure 1).The 5 sites covered a wide range of climatic and ecological characteristics and included pine stands with different specific compositions and characteristics, representing the environmental variations of each forest type (Table 1).In Beteta (BE), covered by Pinus nigra Arn.subsp.salzmannii, and El Pozuelo (EP), a mixed forest of Pinus nigra (84 ± 14%) and Pinus pinaster Ait.(14 ± 14%) (both in the Cuenca Mountains, Iberian Mountain Range, Central Spain), 135 trees were selected for studying (n = 9 plots in each site).In Quintana Redonda (QR, Soria, Iberian Mountain Range, Central Spain), 54 trees (n = 5 plots) in a pure Pinus pinaster stand were selected for studying.In Coirego (CR, Pontevedra, NW Spain), 18 trees in a pure Pinus pinaster stand (n = 18 plots) were selected for studying.In Granadilla (GR, Cáceres, western Central Mountain Range), 32 trees in a mixed stand of Pinus pinaster (82 ± 31%) and Pinus pinea L. (18 ± 31%) (n = 3 plots) were selected for studying.Tree-and stand-level data were thus obtained for 374 individual trees of 3 different species (Pinus nigra, Pinus pinaster and Pinus pinea), across a total 44 plots (12 mixed and 32 pure forest stands).

Prescribed Burning Treatment
The prescribed burning treatments, in order to reduce the vertical and horizontal continuity of the fuel, were carried out by the Forest Service of the respective autonomous regions.In the Mediterranean basin, PB is generally conducted in spring and early autumn.In this context, in Beteta and El Pozuelo, PB were carried out in May and November 2016 and in late June 2019, in different plots for each treatment (in an attempt to simulate high-intensity fires); PB was carried out in Quintana Redonda in March 2015 and April 2018, in Coirego in April 1996, and finally, in Granadilla in November 2015.In all of the treatments, the strip head fire ignition technique was applied downhill, with a head wind.This technique is the most widely used in all areas of study, owing to its flexibility to adapt to different environmental scenarios while enabling sufficient control to produce low-medium intensity fires.This ignition pattern favours the rapid advancement of the front and a shorter residence time of the fire in the soil, thus preventing overheating and excessive consumption of organic matter [92].The litter and shrub layer were the main fire-spread vectors.Air temperature (T, • C), air relative humidity (RH, %) and wind speed (WS, m s −1 ) were recorded during burning at a mobile meteorological station located close to the studied plots.The fire's rate of spread (RS, m min −1 ) was calculated as the mean of several measurements of the time taken for the fire front to travel between two locations separated by a known distance [93].Flame height (FH, cm) was measured from the base of the flame to the top of flame using images taken during the PB and visual observation during PB using metric graduated poles (0.10 m scaled).The main characteristics of the PB are shown in Table 2.

Data Collection and Processing
The data were compiled from experimental fires over a long period of time and by research teams who each used different methods but who all had a common goal of characterizing the thermal regime in the stems and potential cambium damage.The following variables were measured before the PB: stand density (D, trees ha −1 ), diameter at breast height (DBH, cm), diameter at 0.60 m from the ground (D 60 , cm), total height (H T , m), height of the first live branch (H 1 , m), percentage of crown (H 1V , percentage of tree height occupied by crown length), bark thickness (BT, cm) at 0.60 cm from the ground [55] and the coefficient of variation of bark thickness (CV).Bark thickness, defined as the distance from the cambium to the outer bark surface, was measured (to the nearest mm) with a bark gauge.The gauge was inserted into the point of resistance of the cambium.For uneven bark, the maximum thickness was measured at the highest point and the minimum thickness at the lowest point of the bark surface [94].The coefficient of variation of bark thickness was measured as the ratio between the standard deviation and the mean value of the samples.This coefficient was used as a proxy for bark roughness [13].
Potential bole damage was assessed two days after the PB was executed by measuring the maximum and minimum heights of the charred stem (S MX and S MN , cm) and the percentage height of stem charred, relative to the total height (S V , %).The maximum height of charred stem was reached on the leeward side of the stem [7,70,95].Crown damage was assessed using three methods.The observed volume of crown scorch (CS V , %) was estimated visually to the nearest 5%, except in clearer situations where it was possible to make estimates to the nearest 1% [57].Crown scorch height (CS H , m) was calculated using an ultrasonic-based device [96].Litterfall (L, kg ha −1 ) was collected over one year after the execution of the PB in the Cuenca Mountain plots (BE and EP plots, Figure 1, Table 1, n = 9).In order to guarantee the quality and quantity of the litterfall samples, the collection system was designed following the recommendations and parameters outlined in the manual published by the United Nations Economic Commission for Europe (UNECE) within the project entitled "International Co-operative Program on Assessment and Monitoring of Air Pollution Effects on Forests" (ICP Forests, Level II plots) (for further details, see [31,91]).The absolute maximum cambium and outer bark surface temperature (TC MX and TB MX , • C), the duration of temperature above 60 • C (tC 60 , s) at the vascular cambium and the duration of temperature above 300 • C at the bark surface (tB 300 , s) were measured during the PB.The variable tC 60 may be a good indicator of the necrosis of cambium and phloem tissue caused by heat-related damage.A temperature threshold of 300 • C has been found in laboratory-and field-based studies to serve as a proxy for the presence of flame combustion [13,97].Two thermocouples were placed on the leeward side of the trunk of each selected tree, one beneath the bark, at the level of the cambium, and the other on the outer surface of the bark, at a height of 0.6 m from the ground [55].At each location, a hole was drilled through the tree to insert the thermocouple at the cambium level.The thermocouples (type K 1 mm diameter Inconel-sheathed, with a response time of 0.3 s) were connected to data loggers, which recorded the temperature every second [13].

Statistical Analysis
The duration of temperature above 60 • C in the cambium (tC 60 ) is the response variable used as a proxy of cambium damage.Three types of predictor variables were considered (Table 3): (i) variables related to tree and stand attributes; (ii) variables related to fire prescription (burn season); and (iii) variables related to fire severity, comprising a group of variables linked to external signs of damage to the crown and potentially to the stem, and another group reflecting the degree of the heating of the outer bark.Exploratory data analysis (EDA) was initially carried out in order to summarize the main characteristics of the data set.EDA was a critical first step in determining what the data might indicate under the starting assumptions.Covariates with high collinearity were analysed independently to select the most appropriate, not only from a mathematical point of view, but also in a biological sense (this happened, for example, with the bark thickness and the diameter at 60 cm from the base).Scatter plots were constructed to map each variable against an x-or y-axis coordinate, and the relationships between potentially explanatory variables were established (Figure 2).We proposed the application of a random effects hurdle model framework to estimate the probability of the temperature remaining above 60 • C.This approach is commonly used for any form of clustered count data with excess zeros (Figure 3).The model is expressed in two parts: the first part is typically a binary logit model, which models whether the time observed above 60 • C is equal to zero; and the second part determines the time during which that temperature threshold was exceeded.In order to allow for dependence between the two parts of the model, we included parameters which, if taking values other than zero, indicate that the two parts are dependent: a simple classical test can be used here [98].
Hurdle models [99,100] have been developed for modelling zero inflation when regular count models, such as Poisson or negative binomial models, are unrealistic.These models have become popular in many fields, including ecology and environmental science [101][102][103][104]. Hurdle models can be viewed as two-component mixture models consisting of a zero-mass component and a positive observations component following a truncated count distribution, such as truncated Poisson or truncated negative binomial (NB) distribution.
The general structure of a hurdle model is given by Equation (1): where Y i denotes the observed response; i = 1, . . ., n, and n denote the total number of observations; pi is the probability of a subject belonging to the zero component; p(y i ; µ i ) represents the probability mass function (PMF) for a regular count distribution with a vector of parameters µ i ; and p(y i = 0; µ i ) is the distribution evaluated at zero.It can be seen that the positive count is governed by a regular count distribution, as the PMF divided by one minus the PMF of this regular count distribution equals zero.In the present case, the count distribution follows a Poisson distribution; the probability distribution for the hurdle Poisson model is expressed as follows (Equation ( 2)): We implemented two-part zero-inflated mixed-effects models in the novel Generalized Linear Mixed Model (GLMM) adaptive package, which fits the GLMM for a single grouping factor under maximum likelihood by approximating the integrals over the random effects with an adaptive Gaussian quadrature rule [105].This approach allowed us to account for the data structure (zero-inflated, right-skewed continuous data, GLMMzi) and assess the relationship between tC 60 and the set of explanatory variables (Table 3).

DBH
Diameter

Model Assessment and Evaluation
Fitted models were compared according to goodness-of-fit measures, which describe the performance of fitted models for a given data set.The best model was selected on the basis of the Akaike Information Criterion (AIC; Equation (3)) and the Bayesian Information Criteria (BIC; Equation ( 4)).
where df represents the degrees of freedom of fit, and n is the total number of observations in the data.Both the AIC and the BIC impose penalties (stronger with BIC) for the addition of variables to the model, as this increases the sampling variability.Low-AIC and -BIC values indicate that the model provides a good fit for the study data [106,107].

Exploratory Data Analysis
The exploratory data analysis showed an inverse relationship between bark thickness and duration of temperature above 60 • C in the cambium (Table 4) during spring PB, with the same tendency during the summer treatment; whereas, in low-intensity autumn burns, the bark thickness is less determinant (Figure 2a).The highest tC 60 values were reached during the summer burning.For the spring and autumn burning, the data usually do not exceed 500 s of exposure to the 60 • C threshold value.Regarding the variables representative of crown damage, such as crown scorch volume and litterfall, the CS V variable is significantly related to the response variable, so that a higher CS V is directly associated with a higher exposure to temperatures above 60 • C in the cambium (Table 4).During the autumn burning, the tree crowns were hardly affected (Table 2), so the corresponding line does not appear in Figure 2b, showing the low weight of AB data for this variable.Similarly, the tree crowns were scarcely affected during spring burning.As expected, the tree crowns were more affected during summer burning (Table 2), showing the direct relationship of the crown damage to the response variable tC 60 (Figure 2b).A post-fire increase in litterfall does not seem to be related to the duration of temperature above 60 • C in the cambium (Table 4, Figure 2d).A higher height of the first live branch is associated with an increased tC 60 in the studied plots (Table 4; Figure 2c).As expected, a higher mean absolute maximum outer bark temperature implies a greater tC 60 (Table 4, Figure 2g).Regarding the variables associated with stem damage (zero hurdle model coefficients; Table 4), higher values of the maximum height of the charred stem meant that the response variable, tC 60 , was close to 0. Therefore, there is no relationship between a high S MX and a high duration of temperatures above 60 • C in the cambium.Conversely, a high percentage height of charred stem implies either a high tC 60 or that the response variable is not close to 0 (Table 4, Figure 2f).Pinus nigra seems to reach lower tC 60 values (Table 4); although, the data obtained in the studied plots for Pinus nigra are more numerous than Pinus pinaster or Pinus pinea.The tC 60 variable reached a wider range of values during the summer burns (Figure 2h).Regarding stand structure, mixed or pure (Table 4), Figure 2i shows a more normalized distribution of the data in the mixed stand, with lower tC 60 values in the pure stand.

Hurdle Model
The final model accounts for excess zeros in the measured response data.Estimated regression coefficients corresponding to the factors for the best performing zero-inflated Poisson (ZIP) mixed model (according to the AIC value) are shown in Table 4.The intercept of the random effects indicates a decrease in the probability of higher tC 60 values being reached during the autumn PB, and the opposite impact of the treatments applied during summer and spring.The Poisson (count) part of the ZIP model exhibits the risk of a greater duration of exposure as H 1 , CS V and TB MX increase.Furthermore, L and BT have opposite effects to the other variables.Table 4 also includes information obtained from the logistic part of the ZIP model.This part of the model shows the probability of the temperature not reaching values higher than 60 • C. Higher values of S MX reduce the probability of cambium exposure to temperature above 60 • C. On the contrary, higher S V values increase the probability of passing that threshold, with the latter having similar values at the intercept.
For an explanation of these abbreviations, see Table 3.

Discussion
The processes of damage to the cambium in tree stems may be particularly important in the context of PB where low-intensity fires do not generally cause large amounts of crown damage [108].Significant levels of cambium and phloem damage, and the subsequent secondary effects, could reduce the ability of some species to recover after the disturbance [109].This, in turn, could decrease the recurrence of PB or other silvicultural treatments (e.g., thinning), with the stand becoming less resilient to potential wildfires [110].In addition, biotic agents, such as bark beetles and other opportunistic pests, will respond to the decreased tree vigour [111][112][113] and may potentially attack the surrounding live trees [57,89].Consequently, an early diagnosis of post-fire tree damage is of great interest for managers to design adequate fire prescriptions, while enhancing ecosystem resilience.This study's findings showed that the probability of a longer duration of temperature above 60 • C in the cambium in surviving trees depends on the burn season (Table 4).This variable explains the variability of the trees' responses and is useful because of its simplicity, being easily controlled in PB plans [31].The findings showed significant differences in tC 60 between the autumn PB (lower values) and the spring and summer PB, as expected, according to the respective fire behaviour (Table 2, Figure 2).These differences can be explained mainly by seasonal phenology or different fire intensity and fuel consumption levels [114].Overall, spring PB is considered more disruptive to trees because carbohydrate reserves are at their lowest levels at the beginning of the annual growth period [115][116][117][118]. Furthermore, a greater level of damage to fine roots (abundant during this period) can be expected [72,119].However, fuels are usually described as typically drier later on in the year [120], i.e., during the typical burning season (autumn), which may imply higher mortality rates [114,121,122].However, other factors such as the activity of insect pests, the climatic conditions following the fire and tree species susceptibility can interact and lead to a complex response [123] and inconclusive results, probably because fire effects often outweigh the effects of other factors [124].
Considering the existence of stable stand ecotones of Pinus nigra-Pinus pinaster in the Cuenca Mountains and of Pinus pinaster-Pinus pinea in the western Central Mountains Range, pure and mixed stands of the three species were selected throughout the climatic gradient.The model's findings suggested that the temperature remains above 60 • C in the cambium for a shorter time in pure stands than in mixed stands.This could be due to the fact that the species mix was imbalanced, with the dominant species present in a proportion of around 80% (84% in El Pozuelo and 82% in Granadilla).Future research should include stands with a more balanced proportion of species, or even stands with species of different genera.Pinus nigra appeared to be more resistant to low-intensity fires, with less time spent above 60 • C than in Pinus pinaster, although the weight of Pinus pinaster and Pinus pinea in the model is too low to be a deciding factor.Despite this, P. nigra and P. pinaster are considered resistant to low-intensity fire [6,125,126].In general, these species have thicker bark (especially in the lowest part of stem) than other Iberian pines [13,125], such as in the studied plots, where bark thickness is between 1.8 and 4.8 cm.Furthermore, they form relatively open stand structures [6,125] and the height of the first live branch is particularly high [125], being between 6.4 and 9.8 m in the studied sites.However, it must also be pointed out that Pinus nigra and Pinus pinaster forests are among those predicted to be seriously affected by climate change within the Central Iberian Peninsula [127][128][129][130].These stands are therefore more likely to be burned, and determining fire resistance in these species could contribute to improving their management.
Stand variables have also been reported to be good descriptors of fire resilience (e.g., [66,67,69,70,131]).These variables are easy to measure; most of them are typically measured in forest inventories for management purposes or during the monitoring of permanent plots for research purposes [12].In this study, a higher height of the first live branch was associated with an increased tC 60 .This model result is probably related to the structure of some of the stands studied (e.g., high stand density above 1200 tree ha −1 in Beteta).Stand density can strongly affect tree stress, i.e., in dense stands (basal area) where trees must compete for water and nutrients [132], the trees have a smaller diameter and thinner bark, and fire disturbance is expected to be greater [114].Nevertheless, the height of the first live branch in the studied species is high (between 6 and 10 m from the base, which may be compatible with an adaptive characteristic of the species to surface fire) and in most cases, the height of charring was lower than the height of the first live branch, so that the impact on tree survival was not expected to be significant.However, the height of the first live branch must be considered when planning PB in stands with a low height of the first live branch, or, in the case of more severe fire, as in the summer burning, in which a low height of the first live branch appears to be related to cambium damage (Figure 2c).
These findings confirm the key role of bark thickness in protecting the cambium against fire, which is consistent with the findings of numerous previous studies based on empirical and physical approaches, and with recent studies carried out in laboratory and the field involving both mature trees (e.g., [13,30,133,134]) and juvenile trees, and saplings [135,136].Although other bark properties, such as specific gravity, moisture content, bark flammability and physical structure, also influence the resistance of the cambium to heat [10,137,138], it seems that they do so to a lesser extent than bark thickness [139,140].Species that live in fire-prone habitats, particularly those maintained by frequent low-intensity surface fires, generally have thick bark [125,141] that can effectively protect them from the direct impact of low-level fire (as happened during the autumn burning), and other variables related to fire behaviour may have a greater influence on cambium damage [13].However, during intense burning (especially during spring), thinner bark can lead to high tC 60 values, with subsequent damage to the cambium.This could probably be extended to summer burns, although scarce data were obtained in the present study, mainly due to the legal restrictions to carry them out.In most conifers, crown damage is considered the most important factor for determining fire severity level and the best predictor of post-fire mortality at tree level [17,81,131,142].The apparent relationship between tC 60 and crown damage indicated by our empirical model is consistent with previous findings, especially notable in more severe burns.Crown damage in conifers is usually assessed as either the scorch height, assuming that all crown tissue (foliage, buds and twigs) within the scorch height zone is dead [143], or as the volume of crown scorched [131,144].However, some species are less susceptible to thermal damage in the crown [88].This applies to Pinus pinaster, in which needle necrosis occurs after long exposure to temperatures between 55 and 65 • C [125] and the large buds, shielded by scales and long needles, may endure very high levels of defoliation [145].Indeed, P. pinaster ecotypes can display marked differences in post-fire survival for the same volume of crown scorched and mortality, with the probability of death increasing sharply with the number of basal stem quadrants in which cambium death occurs [69].In addition, differences in crown condition can take days to months to become fully apparent to observers, and a delay in quantifying the effects of fire is required, especially considering the possible error associated with the visible estimation of crown scorch volume [88].Litterfall biomass could substitute, or at least complement, commonly used crown damage descriptors.Although this assessment does not require any complicated infrastructure, it requires monitoring for a long time after burning, and therefore this variable may not be readily available.While litterfall was a variable initially considered in the proposed model, it did not appear to be related to higher exposures to the threshold temperature of 60 • C in the cambium (Table 4), probably because data were only available for the Beteta and El Pozuelo plots in the Cuenca Mountains (Table 2).In addition, because litterfall is a stand-scale variable and stem heating is a tree-scale variable, the results would also demonstrate that crown and stem damage are independent if they do not occur together in the same tree [69,146].Beyond these findings, remote sensing techniques show promise for improving estimations of crown damage [147,148].If heat-damaged needles are not capable of re-translocating or reabsorbing labile nutrients, carbohydrate stores and hormones, particularly in ecosystems frequently affected by fire, a cascade of biogeochemical fluxes could occur at the stand or ecosystem level [88].
Stem char height and percentage stem char height are often used as simple proxies for potential cambium damage [16,149,150], although they do not always yield reliable results.In the studied plots, the percentage height of stem charred relative to the total height seemed to be more related to higher tC 60 values.The height of charring is usually higher on the leeward side of the stem, in line with higher combustion turbulence and higher heat generation.This frequently leads to deeper damage on the leeward side, although it may not be sufficient to cause tree death if damage does not occur around the entire stem [26].In this regard, Ref. [151] found that Pinus halepensis trees survived even when 80% of the cambium circumference was killed by fire in the absence of any damage to the foliage.In addition, Ref. [146] observed a lack of any impact on radial growth and physiological variables after death of the cambium in 69% of the stem circumference of P. pinaster trees without crown damage.The above findings suggest a strong resistance to stem damage in both species.Furthermore, in a study conducted three years after fires, Ref. [69] did not observe any mortality in P. pinaster trees in which cambium death had occurred in four quadrants of the trunk, provided there was no crown damage caused by the fire.However, this situation changed dramatically in the presence of bark beetle attacks, with mortality frequently associated with death of the basal cambium.These findings encourage the use of PB in adult trees of both species under an appropriate prescription window and while avoiding damage to the crown and roots.Ref. [112] also proposed that conifer forests are killed by a combined effect of more than 50% bole charring and more than 90% crown scorch, whereas other researchers have indicated that almost 100% girdling is required to kill the tree [7].In any case, this study's findings highlight the importance of considering the percentage stem char height with reference to the duration of temperatures above 60 • C in the cambium; managers should take this factor into account in their burning prescriptions, avoiding charring of a high percentage of the stem (to prevent affecting the cambium) and also stem girdling caused by heat.
The proposed model established a significant relationship between the duration of lethal temperatures in the cambium and the maximum bark temperature (Table 4, Figure 2g).However, these results must be considered with caution, because occasional heating of the bark will not necessarily induce a response in the cambium area.In this respect, the temperature at a specific location within a tree will be influenced by the distribution of temperature over the whole bark surface, which is very irregular under fire conditions [55,152].Only local conditions at the tree level (e.g., the natural accumulation of branches fallen around the tree or lichens growing along the stem) can generate burn intensities that will potentially cause this type of damage [121].This has already been observed in previous field studies investigating the effect of PB on cambium [13].

Conclusions
The proposed model confirms the initial hypothesis, which is that the probability of a longer duration of temperatures above 60 • C in the cambium depends on tree and stand attributes, fire severity descriptors, stem external heating variables related to fire behaviour, and burn season.This model selects variables that are relatively easy to measure in the field and are measured in regular inventories, so as to help managers to elaborate more accurate burning prescriptions.
This study revealed differences in the exposure of the cambium to the threshold value of 60 • C between burning carried out in spring and summer and burning carried out in autumn, with less damage caused in the latter case; although, summer burns were only carried out in Beteta and El Pozuelo, mainly due to the common legal restrictions to carry them out In the case of species with thin bark, which are therefore more susceptible to high temperatures being reached in the cambium, intense prescribed burning should be avoided in the spring or late spring.In this respect, Pinus nigra exhibited a higher resistance, with a shorter time during which temperatures remained higher than 60 • C in the cambium, than P. pinaster or P. pinea, probably because a greater number of P. nigra trees were sampled than P. pinaster or P. pinea trees.However, it is generally known that these three species have different traits enabling them to resist low-intensity fires.Cambium damage is dependent on the stand's structure as a function of the height of the first live branch.This study's findings confirmed the influence of bark thickness in preventing high cambium temperatures during the spread of fire.Crown scorch volume is an important factor for determining the damage to the cambium, although care must be taken during visual estimation of this variable.Litterfall biomass can complement the commonly used crown damage variables, although this study did not reveal a correlation between the litterfall biomass and higher exposures to the threshold temperature of 60 • C in cambium, probably because data were only available for two of the studied areas.
Management actions should aim to avoid high percentages of stem char height, to reduce the potential impact on the cambium, and also stem circumference girdling caused by fire.In addition, management actions should prevent the accumulation of large amounts of fuel around trees as this can lead to high local temperatures of the trees, which have a direct effect on the cambium.

Figure 2 .
Figure 2. Relationships between the main variables selected by the model (Table3) and the response variable (duration of temperature above 60 °C in the cambiu tC60).BT: Bark thickness; CSV: crown scorch volume; H1: height at which first live branch appears; L: total litterfall one year after burn; SMX: maximum height charred stem; SV: percentage height of charred stem; TCMX: mean absolute maximum cambium temperature; AB: autumn burning; MB: summer burning; S spring burning.

Figure 2 .
Figure 2. Relationships between the main variables selected by the model (Table3) and the response variable (duration of temperature above 60 • C in the cambium, tC 60 ).BT: Bark thickness; CS V : crown scorch volume; H 1 : height at which first live branch appears; L: total litterfall one year after burn; S MX : maximum height of charred stem; S V : percentage height of charred stem; TC MX : mean absolute maximum cambium temperature; AB: autumn burning; MB: summer burning; SB: spring burning.

Figure 3 .
Figure 3. Distribution of the tC60 variable (a) with the presence of zeros; and (b) without the presence of zeros.Figure 3. Distribution of the tC 60 variable (a) with the presence of zeros; and (b) without the presence of zeros.

Figure 3 .
Figure 3. Distribution of the tC60 variable (a) with the presence of zeros; and (b) without the presence of zeros.Figure 3. Distribution of the tC 60 variable (a) with the presence of zeros; and (b) without the presence of zeros.

Table 1 .
Main characteristics of the study sites according to species (mean and standard deviation).
BE:Beteta; EP: El Pozuelo; QR: Quintana Redonda; CR: Coirego; GR: Granadilla; D: plot density; DBH: diameter at breast height; D 60 : diameter at 0.6 m from the base; H T : total height; H 1 : height of the first live branch; H 1V : percentage of tree height occupied by crown length; BT: bark thickness; CV: coefficient of variation of bark thickness; "--": data not available.The method of calculating the variables will be explained in Section 2.3.

Table 2 .
Main parameters measured during and after prescribed burning in the study sites (mean and standard deviation).TC MX : mean absolute maximum temperature of the cambium; TB MX : mean absolute maximum temperature of the bark; tC 60 : duration of temperature above 60 • C in the cambium; tC 60C : duration of temperature above 60 • C in the cambium, considering only the trees in which this temperature was exceeded; tB 300 : duration of temperature above 300 • C at the bark surface; tB 300C : duration of temperature above 300 T: air temperature; RH: air relative humidity; WS: wind speed; RS: rate of spread; FH: flame height; S MX : maximum height of charred stem; S MN : minimum height of charred stem; S V : percentage height of charred stem; CS V : crown scorch volume; CS H : crown scorch height; L: total litterfall one year after burn; • C at the bark surface, considering only the trees in which this temperature was exceeded; PC 60 : percentage of trees in which temperature was higher than 60 • C in the cambium; PB 300 : percentage of trees in which temperature was higher than 300 • C at the bark surface; BE: Beteta; EP: El Pozuelo; QR: Quintana Redonda; CR: Coirego; GR: Granadilla; "--": data not available.

Table 3 .
Potential explanatory variables considered in the models and (*) variables selected via the best fit model.

Table 4 .
Summary of the statistical results of the best fit model.