Use of Bayesian Modeling to Determine the Effects of Meteorological Conditions, Prescribed Burn Season, and Tree Characteristics on Litterfall of Pinus nigra and Pinus pinaster Stands

Research Highlights: Litterfall biomass after prescribed burning (PB) is significantly influenced by meteorological variables, stand characteristics, and the fire prescription. Some of the fire-adaptive traits of the species under study (Pinus nigra and Pinus pinaster) mitigate the effects of PB on litterfall biomass. The Bayesian approach, tested here for the first time, was shown to be useful for analyzing the complex combination of variables influencing the effect of PB on litterfall. Background and Objectives: The aims of the study focused on explaining the influence of meteorological conditions after PB on litterfall biomass, to explore the potential influence of stand characteristic and tree traits that influence fire protection, and to assess the influence of fire prescription and fire behavior. Materials and Methods: An experimental factorial design including three treatments (control, spring, and autumn burning), each with three replicates, was established at two experimental sites (N = 18; 50 × 50 m2 plots). The methodology of the International Co-operative Program on Assessment and Monitoring of Air Pollution Effects on Forests (ICP forests) was applied and a Bayesian approach was used to construct a generalized linear mixed model. Results: Litterfall was mainly affected by the meteorological variables and also by the type of stand and the treatment. The effects of minimum bark thickness and the height of the first live branch were random. The maximum scorch height was not high enough to affect the litterfall. Time during which the temperature exceeded 60 ◦C (cambium and bark) did not have an important effect. Conclusions: Our findings demonstrated that meteorological conditions were the most significant variables affecting litterfall biomass, with snowy and stormy days having important effects. Significant effects of stand characteristics (mixed and pure stand) and fire prescription regime (spring and autumn PB) were shown. The trees were completely protected by a combination of low-intensity PB and fire-adaptive tree traits, which prevent direct and indirect effects on litterfall. Identification of important variables can help to improve PB and reduce the vulnerability of stands managed by this method.


Introduction
The projected changes in climate are likely to affect fire regime and increase the fire risk in some Mediterranean countries (including Spain) [1]. Fuel reduction treatments such as prescribed burning (PB) reduce surface continuity and ladder fuels, increase the height of live crown, decrease crown density, retain large trees of fire-resistant species [2], improve success of fire control [3], and may even be effective in creating fire-resilient stands. Nevertheless, fire often has complex effects on different parts of individual trees, microsite conditions, microclimate, and regeneration [4]. In this study, special attention is given to the effect of PB on tree crowns through the subsequent litterfall.
The relationship between litterfall and decomposition determines the depth of the forest floor layer, which acts as a nutrient reservoir in forest ecosystems [5] and plays a key role in protecting the soil against erosion and improving water infiltration [6]. Litter is also an important source of carbon (C) input to the soil [7]. In addition, litterfall biomass has been proposed as a good indicator of stand productivity [8] and it may provide information about the effects of climate change on forests [9]. Litterfall patterns undergo inter-and intra-annual variation (e.g., [10,11]) including spatial variation. Overall, the amount and quality of litterfall depend on climate variables, geographic factors, stand and tree individual characteristics, disturbance of forest management practices, and to a lesser extent, other external variables.
The relationship between litterfall and the duration of drought conditions has been demonstrated in several studies [6,[12][13][14]. In addition, drought conditions may have stronger impacts on Mediterranean ecosystems where rainfall is irregularly distributed [15]. Overall, litterfall production has been reported to be related to latitude, mainly involving temperature and water availability [11]. Indeed, [16] defined a wide range of 3000-11,000 kg ha −1 year −1 of litterfall in different forest types studied worldwide (over a very wide latitudinal range).
Much of the variation in litterfall response to burning can be attributed to stand characteristics (age, structure, species, etc.) [12,17] and to the variable heat sensitivity of different tissues and species [18]. Pinus nigra Arn. ssp. salzmannii (Spanish black pine) and Pinus pinaster Ait. (maritime pine), the species considered in this study, has evolutionary adaptations that help them to persist in fire-prone environments [19,20]. Thus, Pinus nigra has a thick bark, a high crown base height [21], and a self-pruning strategy [22,23]. Bark is considered as the main protective tissue of the cambium [24] and is necessary for photosynthate transport to the crown and nutrient and water storage [25]. In addition, low-intensity PB will burn the less productive part of the crown in larger trees, but scarcely affect the photosynthetic production in the upper part of the crown [26,27]. Maritime pine also has thick bark and high crown base height. In addition, the large buds (shielded by scales and by long needles) may enable very high levels of defoliation, and the storage of seed in serotinous cones facilitates reproduction after fire [19]. In the study area, Spanish black pine often appears mixed with maritime pines. Species diversity has been suggested to promote ecosystem resilience in response to changes through diversity [28].
The use of prescribed fire can influence the patterns and regimes of the litterfall biomass [29]. Indeed, the immediate impact of PB on litterfall and the subsequent recovery of fuel load and structure are considered critical for evaluating the treatment [30,31]. To apply PB safely and minimize the impacts on the ecosystem, weather and fuel related variables must be within the established threshold limits. The burning season is the most controllable feature of the PB regime from a management point of view. However, few studies have investigated the effect of burning season on litterfall. Low-intensity PB may have scarce effects on tree crowns. However, if the PB is very intense, it may cause an increase in litterfall, thereby increasing the risk of a future fire.
The variable scorch height is frequently measured in the field and may be associated with total or partial cambial death or stress, caused by temperatures above the critical threshold [32] or long flame residence times [33]. This could damage phloem and xylem tissues, thus disrupting translocation of photosynthates to the roots [34]. This effect becomes more important as tree height increases because living tissues are less protected by the bark than they are at the base [35].
In a previous study [29], we evaluated the short-term effect of PB in a mixed stand of Pinus nigra and Pinus pinaster and in a pure stand of Pinus nigra. The main aim of the PB conducted was to reduce the accumulated fuel far beyond 50%, which was achieved. The findings showed an increase in the amount of litterfall biomass immediately after PB in spring, especially in the pure stand. The impact of the disturbance was clearly mitigated in the mixed stand. The effect of the treatment decreased over time. In a second study [36], the effects on litterfall of autumn prescribed burning were integrated and results of all treatments in the medium term were analyzed; in addition, information on the nutrient content of litterfall was included. In this regard, the effects (increment in the amount of litterfall) of PB in autumn were scarcely noticeable in the mixed stand and the effects were delayed in the pure stand. The complexity of the litterfall process requires the application of mathematical models to simplify and understand it. The main approach used in the aforementioned study was frequentist inference. However, some of the assumptions of p-values and maximum likelihood estimators are not appropriate for the complex functioning of the litterfall process. A new approach through a Bayesian statistic could solve this problem because Bayesian analysis does not rely solely on a point estimate of the parameter, but rather an entire distribution of possible values. In addition, the Bayesian model may perform better for small sample sizes (also, Bayesian approaches can incorporate more parameters than data) [37]. In this study, this is a key point as litterfall dynamics are strongly dependent on geoclimatic and species characteristics and usually the sample size is small to apply a frequentist approach.
From the prior knowledge of post PB litterfall dynamics, we hypothesized that litterfall may be affected by three broad types of variables: meteorological conditions; stand characteristics (mixed and pure stand) and tree traits that influence fire protection; and fire prescription and behavior. The Bayesian approach involved the analysis of the prior information of [29], data from the other studies cited (e.g., [6,12]) and consideration of the accumulated empirical experience. The aims of the present study were as follows: (i) to ascertain the influence of meteorological conditions after PB on litterfall biomass and to explore the most important variables describing the process; (ii) to explore the potential influence of stand characteristic (mixed and pure stands) and tree traits that influence fire protection (bark thickness and crown base height); and (iii) to assess the influence of fire prescription (season) on litterfall biomass dynamics and fire behavior through easily measurable field variables.

Study Site
Two different sites in the northwest of the Cuenca Mountains (Iberian System) were selected for this study: a mixed stand (MS) of Pinus nigra and P. pinaster in El Pozuelo and a pure stand (PS) of P. nigra in Beteta. According to the data for the last 21 years recorded at the nearest weather station (Cañizares: 940 m a.s.l., provided by the Spanish Meteorological Agency (AEMET) [38]), the mean annual temperature was 12.1 • C (13.2 • C for the period 2016 to 2018) and the average precipitation was 717 mm (599 mm for the period 2016 to 2018). The main characteristics of the two sites were described in [29,36] and a brief summary is shown in Table 1.

Experimental Design
A total of nine square plots of 50 m × 50 m were chosen at each site by applying a randomized complete block design (N = 18). For data collection, subplots of 30 m × 30 m were established to prevent the edge effect ( Figure 1). A total of three treatments were established per experimental site (NB: no burning, SB: spring burning, and AB: autumn burning), with three replicates per block.
All trees in the plots were identified, and the following parameters were measured in each tree: total height (Ht, m), height to the first live branch (H1lb, m), diameter at heights of 0.6 and 1.3 m from the base (D60 and D130, cm), and maximum and minimum bark thickness at 0.6 m from the base (THMxB60, THMnB60 cm). The main stand characteristics of no burning, spring burning, and autumn burning plots in mixed and pure stands are shown in Table 2. (*) Data from [39], (**) According to the classification stablished by [40].

Data Analysis
The Bayesian methodology specifies a probability model that applies prior knowledge about a research parameter, consequently, it is conditioning to the model to perform the adjustment of assumptions. In ecology and environmental analysis, the response variable is commonly measured more than once for each subject across levels of one or more factors, referred to as repeated measures. These types of data cannot be analyzed by linear regression because the residuals for each variable  Note: (*) Data from [29], (**) Data from [36], S, type of stand; PT, plot treatment; Dt, total density; Pn, percentage of Pinus nigra; Pp, percentage of Pinus pinaster; Ht, total height; H1lb, height at which first live branch appears; D60, diameter at 60 cm from the base; D130, diameter at 130 cm from the base; THMxB60, maximum bark thickness at 60 cm from the base; THMnB60, minimum bark thickness at 60 cm from the base; G, basal area.

Prescribed Burning
Prescribed burning was conducted by the Castilla-La Mancha Regional Forest Service in order to reduce the vertical and horizontal continuity of fuel. The spring PB was carried out in May 2016, while the autumn PB was carried out in November 2016. The strip ignition technique was applied at a distance of 1-2 m downhill and with a head wind. This ignition pattern favors rapid advancement of the front and a shorter residence time of the fire in the soil, thus preventing overheating, excessive consumption of organic matter, and high temperatures being reached [41]. This technique is the most widely used in the study area to produce low-medium intensity fire. The efficacy of the burnings was moderate-high (reduction of 59-77% of litter biomass). The main goal of the treatment was therefore successfully achieved [29]. Precipitation (Ortrat, S.L.; KW 3-02), wind speed (Casella; 178031C-3), temperature, and relative humidity (Geonica; STH-5031) were recorded every 10 min at a meteorological station adjacent to the study plots. During the burning, the temperature of the bark (air-surface) and the cambial region (inner bark) of 15 randomly selected trees per plot was monitored at a height of 0.6 m with type K 1 mm diameter inconel-sheathed thermocouples (0.3 s of response time). The thermocouples were connected to dataloggers (DT-USB TCDirect ® ), which recorded the data every second. Variables concerning temperature reached and residence time of flame in bark surface and cambium area were calculated. The threshold value of 60 • C in the cambium area corresponds to the commonly accepted lethal temperature for tree cells [32], although long duration of temperatures of 40-50 • C may cause stress in trees [24]. Maximum and minimum scorch heights were measured after burning. Data collected during and after PB are shown in Table 3. Table 3. Variables measured during and after prescribed burning in mixed and pure stands (mean and standard deviation). TmMxC, mean maximum cambium temperature; TMxB, absolute maximum bark temperature; TMxC, absolute maximum cambium temperature; TRMxB40, percentage of trees in which the maximum temperature in the bark surface (air temperature) was higher than 40 • C; TRMxC40, percentage of trees in which the maximum temperature in the cambium was higher than 40 • C; TRMxB60, percentage of trees in which the maximum temperature in the bark surface (air temperature) was higher than 60 • C; TRMxC60, percentage of trees in which the maximum temperature in the cambium was higher than 60 • C; tB40C, mean time during which temperature is higher than 40 • C in bark; tC40C, mean time during which temperature is higher than 40 • C in cambium; tB60C, mean time during which temperature is higher than 60 • C in bark; tC60C, mean time during which temperature is higher than 60 • C in cambium. For calculating tB40B, tC40C, tB60B, and tC60C, only the trees that exceeded 40 and 60 • C, respectively, were selected.

Litterfall Biomass
Litterfall was collected in fiber glass collectors installed in each plot (eight per plot) after the PB (Figure 1). In order to guarantee the quality and quantity of the samples obtained, the litterfall collection system was designed following the recommendations and parameters outlined in the manual published by the United Nations Economic Commission for Europe (UNECE) under the project entitled "International Co-operative Program on Assessment and Monitoring of Air Pollution Effects on Forests" (ICP Forests, Level II plots) [42]. The spatial distribution of the design covered the entire working area in the 30 m × 30 m plot to ensure the representativeness of the sampling. The specific characteristics of the collectors are detailed in [29]. The litterfall was collected monthly to prevent decomposition of the materials or the chemical leachate. All samples were transported to the laboratory on the same day that they were collected and were oven-dried at 65 • C to constant weight, before being separated into different fractions: needles, cones, inflorescences, branches, bark, and miscellaneous material. The study was carried out between May 2016 and April 2018 in SB plots and between November 2016 and October 2018 in AB plots.

Data Analysis
The Bayesian methodology specifies a probability model that applies prior knowledge about a research parameter, consequently, it is conditioning to the model to perform the adjustment of assumptions. In ecology and environmental analysis, the response variable is commonly measured more than once for each subject across levels of one or more factors, referred to as repeated measures.
These types of data cannot be analyzed by linear regression because the residuals for each variable are correlated. The linear mixed effects model (LMM) is commonly used to analyze repeated measures data. This model combines both fixed and random effects on a linear scale. As the subjects are expected to vary independently in the random effect, correlations between observations within a subject are allowed. Generalized linear mixed models (GLMM) are popular because of their ability to directly recognize multiple levels of dependency and model different types of data. Likelihood-based inference may be unreliable, particularly for small sample sizes, with variance components being particularly difficult to estimate.
In the present study, a generalized linear mixed model approach was used with the aim of establishing the relationship between litterfall and the different variables considered.
The mixed model included 18 plots, each with 24 successive litterfall seizures (µ jk ). The month (k) indicated the time since the PB was carried out (adding the month of the year would only be an unnecessary autocorrelation). There was a lag in the prescribed burnings (spring-autumn) as they were not conducted in the same month. Considering the meteorological variables enabled indirect inclusion of the phenological component. Furthermore, the model considers the order of data collection (i.e., it is a "repeated measures models" application).
We assumed the seizures followed a conditionally independent Gaussian likelihood function: Thus, for the log(µ ij ), a linear model was specified by including the different covariates, x j where α m indicates a set of "fixed" effects for the relevant covariates; u j is a "random" stand effect; and v j is a "random" effect treatment. The statistical inference was carried out by integrated nested Laplace approximation (INLA) implemented in the R-INLA package [43]. INLA uses Bayesian inference on latent Gaussian models by combining Laplace approximations and numerical integration in a very effective manner [44].
The input variables (Table 4) were divided into three groups: variables related to meteorological conditions; stand characteristics (mixed and pure stand) and traits that influence tree protection against fire damage; and fire prescription and behavior. The climatic variables considered (TmMx, TmMm, and Pm) were obtained one month before the litterfall samples were collected [6]. The effect of these variables is not expected to be immediate and the stand response time may be delayed. Cross-validation is commonly used to estimate the out-of-sample prediction error [45,46]. However, researchers require alternative measures such as cross-validation, which involves repeated model fitting and can run into trouble when data are scarce [47]. The index most commonly used for purposes of comparison (e.g., to evaluate which model best fits the data) is the deviance information criterion (DIC) [48,49]. Similarly to the Akaike information criterion (AIC), the DIC has two components: a term that measures goodness of fit, and a penalty term for increasing model complexity. The Watanabe-Akaike information criterion (WAIC) [50] has recently been suggested as a suitable alternative for estimating the out-of-sample expectation in a fully Bayesian approach. This method begins with the computed log pointwise posterior predictive density and then adds a correction for the effective number of parameters to adjust for overfitting [47]. The Watanabe-Akaike information criterion works on the predictive probability density of detected variables rather than on the model parameters and it can therefore be applied in statistical models with non-identifiable parameterization [51]. We also used the conditional predictive ordinate (CPO) [52], which is based on leave-one-out cross-validation to evaluate model assessment. CPO estimates the probability of detecting a value after the others have already been observed. The mean logarithmic score (LCPO) was calculated as a measure of the predictive quality of the model [53,54]. High LCPO values suggest possible outliers, high-leverage, and influential observations. Finally, we used an area under operating curve score (AUC) approach to calculate the predictive accuracy of each method by comparing the validation data with the predicted presence value. The AUC score, a commonly used and adequately performing measure of predictive accuracy [55], calculates the relative numbers of correctly and incorrectly identified predictions across all possible classification threshold values of the binomial response. An AUC value equal to or below 0.5 indicates a predictive ability equal to random expectation, and a value of one indicates a perfect predictive ability [56].

Litterfall
Mean total litterfall production for all treatments was 3054 ± 339 kg ha −1 year −1 in the MS and 3050 ± 674 kg ha −1 year −1 in the PS (Table 5). Intra-and inter-annual variability was observed in both stands (Figure 2a-c). Maximum amounts of litterfall were collected between June and September, representing 63% of the total biomass in the MS and 51% in the PS (mean values for all treatments). Maximum amounts were mainly collected in August in both areas, while minimum amounts were collected in winter months (December, January, and February). Needles were the most important fraction, representing 58% of the total biomass in the MS and 57% in the PS (mean values of all treatments) [36]. There was a short-term effect (3-4 months) on litterfall after SB in the MS and the PS, although the differences were greater in the PS (Table 5). In the following months, no differences were observed in the MS (the litterfall biomass was actually greater in the NB plots). In the PS, the differences remained throughout the study period, although they decreased steadily over time. No effect was observed after autumn-prescribed burning in the MS, but the effect of autumn burning was detected in the PS.

Bayesian Model
Graphical representation of INLA estimation for the fixed effects is presented in Figure 3. This chart presents nine of the 17 variables finally selected by the model (Table 4) and how they are related to litterfall. Variables distributed on the positive side are positively related to litterfall and those on the negative side are negatively related to litterfall. The variables that appear in both positive and negative areas have random effects on litterfall. This distribution is summarized in more detail in Table 6, with mean, standard deviation, and 95% credibility interval for the different variables. The mean values (disks) and 95% credible intervals (lines) highlight the estimated effect of the coefficients of variation. Credible intervals crossing zero (dashed line) suggest that the absence of any effect of the corresponding covariate cannot be ruled out, given the data and the model assumptions. The distribution of random effects is also shown in Table 6. Random effects are often used to account for over dispersion in Poisson models, and high values of random effects show that covariate explains most of the over dispersion in the data. However, if the scale of the estimated random effects is similar to that of the estimated covariates, both the covariates and the random effects will explain the overdispersion.
Meteorological variables (mainly mean maximum temperature, TmMx; days with snow, SwD; days with storm, StD) were the most important variables influencing litterfall biomass. The proposed model provides a positive response concerning the variability in litterfall due to two types of stands. The findings showed that minimum bark thickness and the height of the first live branch had random effects on the response variable. Although the treatment also had a positive effect on litterfall, it had a slightly weaker effect than the type of stand. Scorch height on the leeward side did not affect litterfall.
The variable time during which the temperature was higher than 60 • C in bark and cambium were selected by the proposed model, although the effects of these variables were random.

Discussion
Although the importance of litterfall in the dynamics of the recovery of forest soils is recognized [57], scarce information about the effects of PB treatments on litterfall has been published. When planning prescribed burns, forest managers must reach decisions in the face of uncertain conditions, as litterfall patterns depend on several variables that often interact in a complex way. The Bayesian approach seems to be a good statistical model for dealing with complex models (it enables the inclusion of multiple variables) and for implementing prior knowledge [29,36], thus enabling more robust conclusions to be obtained.
An immediate impact on litterfall was detected after spring and autumn PB. The effect was stronger in the pure stand. Over time, the effect of the treatment decreases. The different effects of PB on the two types of stand strengthen the results obtained in the statistical analysis where the random variable "stand" directly influenced the litterfall. A gradual decline in the effect of the disturbance has also been described by [6,58] in different thinning regimes. Despite this immediate effect, mean total litterfall collected was within the range of natural litterfall in Spanish stands (e.g., [6,12]) and close to those obtained (3337 ± 841 kg ha −1 year −1 ) by [59] in a Pinus nigra stand in Mora de Rubielos (Teruel).
Meteorological variables (in general) were the most important variables influencing litterfall biomass. Indeed, the marked inter-annual variability in climate in most Mediterranean areas drives the inter-and intra-annual variability in litterfall [10][11][12]. Litterfall dynamics may be more sensitive and increase when low-intensity prescribed burning is followed by snow (not significant variable) or storm events. Snowfall and storms have been suggested to be the cause of peaks in the abundance of some fractions of litter such as green needles, branches, and bark (e.g., [12,29,60]). Data provided by AEMET (1997-2018) [38] showed that a range of months with the most snowy days occurred mainly Figure 3. Estimation of fixed effects. TmMx, mean maximum temperature; SwD, days with snow; StD, days with storm; Pm, monthly precipitation; THMnB60, minimum bark thickness at 60 cm from the base; H1lb, height at which first live branch appears; HSmMx, mean maximum scorch height; tC60C, time where temperature is higher than 60 • C in cambium; tB60C, time where temperature is higher than 60 • C in bark. (*) Indicates significant effects. Note: TmMx, mean maximum temperature; SwD, days with snow; StD, days with storm; Pm, monthly precipitation; THMnB60, minimum bark thickness at 60 cm from the base; H1lb, height at which first live branch appears; HSmMx, mean maximum scorch height; tC60C, time where temperature is higher than 60 • C in cambium; tB60C, time where temperature is higher than 60 • C in bark. Values shown in bold indicate significant effects.

Discussion
Although the importance of litterfall in the dynamics of the recovery of forest soils is recognized [57], scarce information about the effects of PB treatments on litterfall has been published. When planning prescribed burns, forest managers must reach decisions in the face of uncertain conditions, as litterfall patterns depend on several variables that often interact in a complex way. The Bayesian approach seems to be a good statistical model for dealing with complex models (it enables the inclusion of multiple variables) and for implementing prior knowledge [29,36], thus enabling more robust conclusions to be obtained.
An immediate impact on litterfall was detected after spring and autumn PB. The effect was stronger in the pure stand. Over time, the effect of the treatment decreases. The different effects of PB on the two types of stand strengthen the results obtained in the statistical analysis where the random variable "stand" directly influenced the litterfall. A gradual decline in the effect of the disturbance has also been described by [6,58] in different thinning regimes. Despite this immediate effect, mean total litterfall collected was within the range of natural litterfall in Spanish stands (e.g., [6,12]) and close to those obtained (3337 ± 841 kg ha −1 year −1 ) by [59] in a Pinus nigra stand in Mora de Rubielos (Teruel).
Meteorological variables (in general) were the most important variables influencing litterfall biomass. Indeed, the marked inter-annual variability in climate in most Mediterranean areas drives the inter-and intra-annual variability in litterfall [10][11][12]. Litterfall dynamics may be more sensitive and increase when low-intensity prescribed burning is followed by snow (not significant variable) or storm events. Snowfall and storms have been suggested to be the cause of peaks in the abundance of some fractions of litter such as green needles, branches, and bark (e.g., [12,29,60]). Data provided by AEMET (1997-2018) [38] showed that a range of months with the most snowy days occurred mainly between December and March. During this period, for all treatments, branches represented on average 40% and 52% (MS and PS, respectively) of the total litterfall collected [36]. Branches may break as a result of the weight of accumulated snow, especially branches that are partly affected by fire. In addition, radiation or convection processes may increase thermal pruning. However, the low-intensity burning carried out in this study does not allow us to confirm the existence of this effect. The months with the most stormy days occurred between May and August. Nonetheless, findings on the influence of storms in litterfall biomass may interact with the time of natural abscission of needles [29].
Prescribed burning must be planned carefully, especially if carried out after months of high temperatures and consecutive months of low rainfall. The mean maximum temperature is a significant variable affecting litterfall in the proposed model. However, the model showed that the effects of mean precipitation varied at random, in contrast to previous findings on the influence of physiological drought on litterfall [6,17]. This may be because the precipitation considered corresponded to the month prior to collecting the litterfall, and its influence on litterfall may be due to periods of extended drought when the accumulated precipitation is below the tree requirements. It is widely known that temperature and water regulate biological processes. Indeed, litterfall biomass is most abundant in summer (June-September) (Figure 2) in the stands under study, representing mean percentages of 63% and 52% of the annual biomass in the MS and in the PS, respectively. Likewise, the combination of increased temperature and longer duration of light period in warmer months increases the activity of growth hormones that cause the abscission of old needles [61]. Information about the ecosystem response to warming is essential for determining the impacts of perturbation in forest systems. An increase in temperature of 1.5 • C is predicted to occur between 2030 and 2052, if temperatures continue to increase at the current rate [62]. The effect of meteorological conditions must therefore be taken into account in considering the impact of PB on litterfall. Other variables indirectly related to meteorological conditions such as phenology may influence the litterfall dynamics, although to a lesser extent [63].
The proposed model provides a positive response concerning the variability in litterfall due to type of stand (mixed or pure). As mentioned, a slight damping of PB effects in the mixed stand was already observed in the previous study [29]; the present results, which cover a longer period of time, reinforce these findings. In the Cuenca Mountains, the Pinus nigra-Pinus pinaster ecotone generates stable stands, which should be taken into account in PB. Exploration of the resilience and vulnerability of mixed and pure stands to perturbations such as PB is essential to establish recommendations for management and fire prevention strategies.
Minimum bark thickness had a random effect in litterfall. In the present study, the PB was of low-intensity (the typical type of burning conducted in these Pinus nigra stands). The bark thickness was therefore probably sufficiently thick to protect trees against thermal damage (even the minimum values measured) [64], thus guaranteeing supply of water and nutrients to the needles (and preventing premature fall). The thick bark of Pinus nigra is considered an adaptation to surface fire [21,65], and this species can persist after being affected by surface fires during several centuries [20]. The same applies to maritime pine populations [66]. As bark thickness is an easily measurable variable and is correlated with the time required to reach lethal temperatures in the ca7mbium [67,68], it can help to predict the resistance of bark to a surface fire and to prevent the crown being indirectly affected.
Although high crown base height may enable trees to escape the effects of surface fires such as those induced by PB, the study findings showed a random relationship between the amount of litterfall biomass and the height of the first live branch. In MS, the first live branch appeared at 6.4 ± 1.8 m from the base, and in PS, at 8.2 ± 2.5 m (Table 1). In most cases, scorch height was lower than the height of the first living branch, which may explain the observed responses in Bayesian models (for low crown base height or more variable values, the relationship may be stronger). It is widely accepted that Pinus nigra exhibits a self-pruning strategy [20,23] with high crown base heights [21] that enable adult trees to resist fire. The present data are consistent with the previously reported crown base height of 7.3 m for the same species in stands in northeastern Spain [22]. Nevertheless, an average base height of 4.0 ± 0.15 m has been reported for Pinus nigra stands in eastern Spain [20]. A low height of the first live branch must be taken into account in planning prescribed burning as it will likely affect litterfall biomass.
Along with type of stand, the proposed model provided a positive response concerning the variability of litterfall due to treatment (NB, SB, and AB). However, type of stand had a slightly stronger effect on litterfall. This may be due to the low intensity of prescribed burning and because the immediate effect of PB (3-4 months after PB) decreased over time (24 months), thus the accumulated litterfall may thus disguise these immediate effects produced by burning. However, if PB has a severe and immediate impact, it may cause an increase in accumulated fallen biomass, thus leading to a higher risk of future fire and a short-lived treatment efficacy [30].
Scorch height on the leeward side had a random effect on litterfall. The maximum scorch height was reached on the leeward side of the trunk (63 ± 70 cm in MS and 107 ± 110 cm in PS, Table 3), although below the mean height of the first live branch. Differences in scorch height on the leeward and windward sides (12 ± 20 cm in MS and 21 ± 40 cm in PS, Table 3) were due to the so-called "chimney effect" caused by the incidence of tree on the flame geometry [69]. Consequently, there are differences between the expected (prescribed) and real height of the flame at which the stem will be heated. Overall, a high scorch is caused by the effect of flame height, which could translate into stronger effects of radiation and convection processes negatively affecting the crown. Nevertheless, scorch marks on the highest parts of the stems in the stands under study were mainly due to the presence of lichens, which burned briefly. Furthermore, no trees were completely scorched, except seedlings under canopy and small dominated trees. In this respect, analysis of the nutrient contents of needle fraction collected reinforced the findings, showing nitrogen (N) concentrations corresponding to old needles in abscission time (4.5 and 5.1 mg g −1 mean value for all treatments and periods in El Pozuelo and Beteta, respectively) [36].
The time during which the temperature was higher than 60 • C in bark and cambium had random effects on litterfall. In PB, in both spring and autumn, the temperature in the cambium area only exceeded 60 • C in 17% of trees in MS and in 8% of trees in PS ( Table 3). As Bayesian models are based on an iterative learning process and the findings depend on previous knowledge about the system and new evidence (from the data), it is reasonable to consider that the study variables will have random effects on litterfall. This may also be explained by the fact that the burning was not completely homogeneous over the whole plot, but depended on fuel accumulation close to the trees. The process of cambial damage may be especially important in the context of PB, in which relatively low intensity fire does not generally result in large amounts of crown damage [68]. However, more data on moderate-high intensity fires are needed to clarify the effects of these variables.

Conclusions
The effects of meteorological variables on litterfall were more significant and even more important than variables associated with stand characteristics or with prescription burning. Mean maximum temperature in the months prior to prescribed burning must be considered in decision making processes, because high temperatures may enhance the accumulation of fallen biomass, decreasing the duration of treatment effects and finally enhancing the fire risk. The final model did not indicate a close relationship between mean precipitation and litterfall. Presumably, an extended period of drought conditions is required to induce stress in trees and increase litterfall.
The variability in litterfall was affected by stand characteristics and tended to dampen the effect of perturbation in the mixed stand. Adaptation of bark thickness and high first live branch height in the species studied was sufficient to provide complete protection against low-intensity fire.
Although the treatment affected the variability in litterfall, the stand characteristic had a slightly stronger effect on litterfall. The results suggested that burning in spring (when the mean temperature is generally higher than in autumn) may have a greater impact on litterfall. In addition, the maximum litterfall production begins immediately afterward and may reinforce the effect. However, one-off episodes of snowfall and storms during winter-autumn may increase the abundance of litterfall, and the effect may even be maintained throughout the year until the time of maximum leaf-fall. Maximum scorch height (leeward size) did not have effects on litterfall. This variable is easily measurable after the passage of fire, and it can help in adjusting prescriptions and in estimating the effects on litterfall. The cambial temperature reached did not affect the live cells, and it therefore had no indirect impact on litterfall. This is the first report of the use of a Bayesian approach to describe the complex litterfall dynamics after PB, with good results. However, the model should be tested with a larger database including information from more moderate-severe fires in order to improve PB planning.