Franck Jabot
email: franck.jabot@irstea.fr

Why preferring parametric forecasting to nonparametric methods?

Keywords: chaotic dynamics, model checking, nonparametric forecasting, observation error, time series analysis

A recent series of papers by Charles T. Perretti and collaborators have shown that nonparametric forecasting methods can outperform parametric methods in noisy nonlinear systems. Such a situation can arise because of two main reasons: the instability of parametric inference procedures in chaotic systems which can lead to biased parameter estimates, and the discrepancy between the real system dynamics and the modeled one, a problem that Perretti and collaborators call "the true model myth". Should ecologists go on using the demanding parametric machinery when trying to forecast the dynamics of complex ecosystems? Or should they rely on the elegant nonparametric approach that appears so promising? It will be here argued that ecological forecasting based on parametric models presents two key comparative advantages over nonparametric approaches. First, the likelihood of parametric forecasting failure can be diagnosed thanks to simple Bayesian model checking procedures. Second, when parametric forecasting is diagnosed to be reliable, forecasting uncertainty can be estimated on virtual data generated with the fitted to data parametric model. In contrast, nonparametric techniques provide forecasts with unknown reliability. This argumentation is illustrated with the simple thetalogistic model that was previously used by Perretti and collaborators to make their point. It should convince ecologists to stick to standard parametric approaches, until methods have been developed to assess the reliability of nonparametric forecasting.

Introduction

The societal demand for predictions in ecology poses a major challenge to the scientific community. It triggers the need to design efficient methodologies that can tackle the inherent complexity of ecosystems. Most predictive exercises in ecology so far have relied on parametric approaches, that is on the design of a model, on the estimation of model parameters that best fit available data, and on the use of the calibrated model to make predictions. Parametric methods encompass widely different approaches to predictions, from correlative approaches, in which the model is a statistical one (e.g., [START_REF] Austin | Spatial prediction of species distribution: an interface between ecological theory and statistical modelling[END_REF], to detailed individual-based modeling approaches [START_REF] Grimm | Individual-based modeling and ecology[END_REF][START_REF] Hartig | Statistical inference for stochastic simulation models -theory and application[END_REF][START_REF] Jabot | EasyABC: performing efficient approximate Bayesian computation sampling schemes using R[END_REF]. In this general methodology, the generally accepted view is that by refining our understanding of the processes at stake in the dynamics of a system, we will improve our ability to make reliable predictions. This is one of the arguments advanced in favor of complex mechanistic modeling in ecology [START_REF] Evans | Do simple models lead to generality in ecology?[END_REF].

In contrast with this general approach, Perretti and colleagues have argued in a series of recent papers that model-free nonparametric approaches may be a better way to tackle ecosystem complexity (Perretti et al. 2013a,b,c). These nonparametric approaches rely on the idea that one should be able to forecast future dynamics by analogy with past dynamics. More precisely, a phase space representation of past dynamics is used, in which variables of interest are plotted as a function of each other, possibly with time delays [START_REF] Kantz | Nonlinear time series analysis[END_REF]. Forecasts are then built using the collection of previously observed dynamics, weighted by their proximity in the phase space (see Materials and Methods and Perretti et al. 2013a,b for additional details). Perretti et al. (2013a,b) showed that such a nonparametric approach to forecasting could outperform parametric forecasts in some stochastic dynamical systems in chaotic regimes, even in situations in which the parametric forecast was based on the same model that was used to generate the data. This poor performance of parametric forecasting was due to two main reasons.

The first reason is the well known instability of parametric inference in chaotic systems [START_REF] Berliner | Likelihood and Bayesian prediction of chaotic systems[END_REF][START_REF] Wood | Statistical inference for noisy nonlinear ecological dynamic systems[END_REF]. [START_REF] Hartig | Does model-free forecasting really outperform the true model?[END_REF] recalled an alternative method to calibrate the model which is based on a segmentation of the time series to circumvent the sensibility to initial conditions encountered in chaotic regimes [START_REF] Pisarenko | Statistical methods of parameter estimation for deterministically chaotic time series[END_REF]. Armed with this trick, [START_REF] Hartig | Does model-free forecasting really outperform the true model?[END_REF] have shown that the parametric forecast is more accurate than the nonparametric one in the case where the same model was used to generate the data and to ground the forecast.

The second reason for parametric forecasting failure was called by Peretti et al. (2013c) the "true model myth", and comes down to the very meaning of what a model is [START_REF] Box | Science and statistics[END_REF]. Indeed, there is a necessary mismatch between reality and its simplified modeled representation. Perretti et al. (2013c) therefore repeated their analysis, using the theta-logistic model (see Methods) to generate a virtual dataset, and basing the parametric forecast on the simpler logistic model. In this setting, and using the inference method proposed by [START_REF] Hartig | Does model-free forecasting really outperform the true model?[END_REF], they found that their nonparametric forecast was more accurate than the logistic based forecast.

As already argued by [START_REF] Hartig | Does model-free forecasting really outperform the true model?[END_REF], parametric approaches offer the advantages of transferability and theoretical understanding which are lacking in nonparametric methods. They further enable to mobilize previous knowledge on the system studied through prior distributions in Bayesian analyses. The aim of this contribution is to demonstrate that parametric approaches are currently also preferable in a forecasting perspective, using the theta-logistic example of Perretti et al. (2013c). Two key comparative advantages of parametric approaches over nonparametric ones will be detailed. First, the likelihood of parametric forecasting failure can be diagnosed thanks to simple Bayesian model checking procedures. Indeed, the convergence of a parametric inference method does not ensure that the fitted model is a satisfactory representation of the data. Simple methods are available to check, with minor computing costs, the ability of the fitted model to recover major data characteristics [START_REF] Rubin | Bayesianly justifiable and relevant frequency calculations for the applied statistician[END_REF][START_REF] Gelman | Posterior predictive assessment of model fitness via realized discrepancies[END_REF]. In cases where the fitted model fails at such model checking tests, modifications of the model should be sought, before any reliable predictions can be made. Second, when parametric forecasting is diagnosed to be reliable, forecasting uncertainty can be estimated thanks to virtual data generated with the fitted to data parametric model. This second advantage enables to complement a forecast with a detailed understanding of its inherent limitations. These two key advantages enable scientists to enrich their forecasts with a description of their reliability.

Methods

The theta-logistic model

In the following, the theta-logistic model with observation error is used to generate time series of population abundances. These time series represent the observed "reality" that we are trying to forecast. The deterministic theta-logistic model is defined by the following difference equation:

X(t+1) = X(t) × r × [ 1 -( X(t) / K ) θ ] (1)
where X(t) is the true population abundance at time t, r is the population growth rate, K is the carrying capacity and θ is a model parameter. When θ equals 1, the model is the well known logistic model of population dynamics, which can generate deterministic chaos [START_REF] May | Biological populations with nonoverlapping generations: stable points, stable cycles, and chaos[END_REF].

Equation (1) describes the true dynamics of the population, while observers only have access to noised estimations of population size N(t). Following Perretti et al. (2013b), the observation error is modeled here using:

N(t) = X(t) × exp[ ε(t) ] (2)
where ε(t) are independent normal random variables with mean 0 and standard deviation s obs , leading to lognormal observation error.

Equations ( 1) and ( 2) are used to generate artificial data. Six values of θ were used {1; 1.1; 1.2; 1.5; 2; 4} with corresponding r values {3.7; 3.7; 3.5; 3; 2.5; 1.8} to ensure a chaotic regime and population survival. Three levels of noise s obs were used: 0.05, 0.1 and 0.2 (represented in Figure S1), leading to an increase of the temporal variance in observed population sizes N(t) of 3, 12 and 52% respectively, compared to the temporal variance of real population sizes X(t). This range of observation noise level is in line with previous theoretical studies ( [START_REF] De Valpine | Fitting population models incorporating process noise and observation error[END_REF] and with reported values estimated on the global population dynamics database [START_REF] Knape | Are patterns of density dependence in the global population dynamics database driven by uncertainty about population abundance?[END_REF]. For each triplet of values (θ; r; s obs ), 50 time series of 100 time steps were generated, starting from an initial population size X(1) equal to 0.5, so that 900 time series were generated in total.

Parametric model fit

Each of these 900 time series was fit with the logistic model, which is the special case when θ equals 1 in equation 1. Only the 50 first steps of the time series were used for calibration, and the 50 remaining steps were used to assess the quality of the forecasts. Following [START_REF] Hartig | Does model-free forecasting really outperform the true model?[END_REF], the time series were split into 10 subsets of 5 time steps each, and the X(t)

values at the beginning of each subset were considered unknown. The rationale of this splitting procedure is to focus on short subset of the time-series so as to decrease the impact of the sensitivity to initial conditions which is specific to chaotic dynamics. The length of the splitting procedure was guided by [START_REF] Hartig | Does model-free forecasting really outperform the true model?[END_REF] who have already shown that it was efficient for this parametric model. For other case studies, alternative splitting length might be more efficient. They could be compared by trial and error. Three MCMC chains were used with dispersed initial r values, 500,000 adaptation steps, 500,000 burn-in iterations and 1,000,000 posterior samples, of which 10,000 (1 every 100 steps) were retained to perform model checks and compute predictions. Prior distributions were uniforms on [2;4.5], [0.01;10] and [0.1;40] respectively for r, K and the observation precision equal to 2 / 1 obs s . This last prior distribution was changed to uniform on [10;600] for cases of low and middle observation noise (s obs equal to 0.05 or 0.1 in the simulations). Convergence of the chains was checked with the multivariate potential scale reduction factor, that was required to be smaller than 1.2. Bayesian analysis was performed with R version 3.0.1 and JAGS 3.4 [START_REF] Plummer | JAGS: a program for analysis of Bayesian graphical models using Gibbs sampling[END_REF]. R and JAGS scripts are provided in the Appendix.

Model checking

The 10,000 parameter sets sampled in the MCMC were used to simulate 10,000 posterior predictive datasets with observation error of 50 time steps (for each of the 900 fitted time series).

Note that this implies to have additionally recorded during the MCMC the fitted initial population size X(1). Five summary statistics were computed on each posterior predictive dataset and on the "real" fitted time series: the average observed population size N , the standard deviation of observed population sizes through time N sˆ, and the autocorrelation coefficients R(t) of observed population sizes with time lag t equal to 1, 2 and 3 (

( ) ( ) ( ) [ ] 2 1 / ) ( ) ( 1 N N t N t N N t N N R s t t t t ∑ + = - - × - - =
). For each summary statistic, a p-value was computed as the percentage of posterior predictive statistics below the value computed on the real dataset. The minimum p min of these 5 p-values (or of 1 -p-value in cases it is lower than the p-value) is reported in the result section.

Posterior predictive p-values are not uniformly distributed between 0 and 1, but instead tend to be peaked around 0.5 [START_REF] Robins | Asymptotic distribution of p values in composite null models[END_REF][START_REF] Gelman | Two simple examples for understanding posterior p-values whose distributions are far from uniform[END_REF] We also tested the method proposed by [START_REF] Crespi | Bayesian model checking for multivariate outcome data[END_REF] to perform model checking on multivariate outputs. It was however found to be less efficient at discarding poor model fits than the above approach, in the present case study. Note also that Peretti et al. (2013c) mentioned that they performed some model-fit diagnostics when applying parametric approaches, but without further explanation. The fact that they were not able to detect model misfit is likely due to the use of too liberal model checking tests. Such too liberal tests are likely to stem from a direct interpretation of posterior predictive p-values or to the use of a set of summary statistics that is not sufficient to properly assess model fit quality. This should serve as a reminder that the power of model diagnostic procedures should be tested with simulated datasets.

Parametric Forecast error

Following [START_REF] Hartig | Does model-free forecasting really outperform the true model?[END_REF], point estimates of the various parameters (the medians of the parameters sampled during the MCMC) were used to compute predicted future population sizes for up to the next 50 time steps. Concretely, for each parameter point estimates, 50 time series were computed starting at time step t (with t in [50;99]), and lasting during 100 -t steps.

These time series were using as initial X(t) values, a maximum-likelihood estimator based on the last five observed population sizes N(t-4)…N(t), and using the model parameter point estimates.

The mean square error between the simulated X(t) values and the reference X(t) values were then Overall forecasting efficiency was summarized in the statistics Area SRMSE defined as the overall improvement of the forecast compared to a forecast equal to the average population size (having a SRMSE equal to 1). Area SRMSE was computed as:

( ) ∑ = - = 8 1 0 ), ( 1 max i SRMSE i SRMSE Area
where SRMSE(i) is the SRMSE of a forecast i time steps in the future. Negative values for 1-SRMSE(i) were not considered in this formula to focus on the short-term forecast efficiency of the parametric and nonparametric methods, since both will fail for long-term predictions in chaotic systems (Perretti et al. 2013b).

Nonparametric forecasting

Nonparametric forecasting was also implemented on these 900 time series using the S-Map model as in Perretti et al. (2013a,b). This consists in first, transforming the univariate time series N(t) into a time series of time-delayed vectors (N(t),N(t-1),...,N(t-E+1)), with E the embedding dimension. If n temporal points are available in the time-series, this produces n-E+1 vectors of dimension E. Each vector (N(t),N(t-1),...,N(t-E+1)) is associated to the following trajectory ((N(t+1),...,N(t+k)), observed in the k following time steps. The nonparametric forecast of the k next step consists in a weighted average of these following trajectories over the n-E+1-k vectors that have k following time steps available in the time series. The weight of a vector is given by ( ) ( )

d d d w / exp θ - =
where θ is a parameter of the S-Map model, d is the Euclidean distance between the current data and the focal vector whose weight is computed, and where d is the average of d over the various vectors. Embedding dimensions E equal to 2, 3 and 4, and θ values between 0.5 and 25 (with an increment of 0.5) were tested. The combination (E,θ) leading to the smallest mean square error on a forecasting interval of length 1 was determined by crossvalidation as in Perretti et al. (2013b). Note again that, contrary to what Perretti et al. (2013a) did, observed population sizes N(t) in the past were used to ground the forecasts, instead of true population sizes X(t), since this is the only piece of information available in real settings. R scripts are provided in the Appendix.

We repeated all the previous analyses including lognormal process errors in the model (with s proc equal to 0.005, and using 10 instead of 50 repetitions per parameter set). Results were qualitatively similar, although MCMC convergence was slower. These results are therefore not reported further in the following.

Results

The model checking procedure allowed to detect, in the majority of cases, the failure of the logistic model (θ=1) to satisfactorily fit the data in cases of low (s obs =0.05) and middle (s obs =0.1) noise level (Figure 1A,B). Indeed, 94% of the p min values obtained on data generated with the theta-logistic (θ =4) were lower than the 5 percentile of p min values obtained on data generated with the logistic model (θ =1), in the low noise case (Figure 1A). This number was still equal to 64% in the middle noise case (Figure 1B). In contrast, this model checking procedure failed in the large noise case (Figure 1C). The ability of the model checking procedure to detect model misfit was also influenced by the discrepancy between the model and the simulated data, as measured by the θ value used to simulate the data (Figure 2). In the low and middle noise cases (Figure 2A,B), model misfit could be detected for values of θ larger than or equal to 1.5, since p min values were sharply smaller than those obtained with the reference model (with θ equal to 1). Model misfit could not be detected in the large noise case, irrespective of the value of θ (Figure 2C). As the discrepancy between the model and the simulated data (measured by the θ value used to simulate the data) rose, forecast errors progressively increased for parametric forecasts, especially in the low noise scenario (Figure 3A, B andC). In contrast, nonparametric forecasts were only marginally affected by the value of θ used in the simulated datasets and were found to outperform parametric forecasts for large value of θ (Figure 3C and F), while the contrary was observed for smaller θ values (Figure 3A, D andE). SRMSE reached larger than one values for both approaches, indicating that they both performed worse than a simple average forecast in the long term. This result is not surprising given the notorious sensitivity to initial conditions inherent to chaotic systems, which prevents any reliable long-term forecast. No reliable predictions could be achieved for more than 5 time steps in the future in all cases. Parametric and nonparametric approaches were comparable, in all but one case (Figure 3C). For this low-noise case with strong model misfit, the nonparametric approach strongly outcompeted the parametric one. Notably, and contrary to what Perretti et al. (2013b,c) argued, in large noise cases, forecast efficiencies were poor for both parametric and nonparametric approaches (Figure 3D, E and F). Finally, overall forecast efficiencies of the parametric approach were found to decrease with the level of observation noise and to be linked with model checking statistics (Figure 4). In both the low and middle noise cases, when model diagnostics would suggest that the model fit quality is acceptable (large p min in Figure 2A, B obtained for θ values lower than 1.5), reasonably good parametric forecast accuracies were reached that were comparable to nonparametric ones (Figure 4A,B). For larger θ values, parametric forecasts were less efficient than nonparametric ones (Figure 4A,B). In the large noise case, forecasting accuracies were poor for both parametric and nonparametric approaches (Figure 4C), and model checking procedures failed at detecting model misfit (Figure 2C). Importantly though, in this last case with large observation noise, forecasting quality was poor even when data were simulated with the model fitted (Figure 4C for a θ value equal to 1). This means that ecologists could in this case detect the lack of parametric forecast reliability, using simulated data with the fitted level of observation noise. In contrast, when forecasts performed with the nonparametric approach fail (Figure 4F), ecologists have currently no way to detect this failure.

Discussion

Why using parametric forecasting?

This study aimed at highlighting two key advantages of parametric approaches to forecasting.

First, Bayesian model checking procedures are powerful tools to asses the misfit of a model calibrated on real data. These tools enable to identify cases in which our modeled understanding of an ecological system is too crude to make reliable predictions. In such cases, modifications of the model should be sought, before any reliable predictions can be made. By re-examining the toy example of the theta-logistic model proposed by Perretti et al. (2013c), it was shown that a simple model checking procedure was sufficient, provided that observation noise was not too large, to distinguish cases in which model misfit was large and parametric forecasting poor from cases in which model misfit was not evidenced and parametric forecasting was reasonably good (Figure 2 A second key advantage of parametric approaches is to provide ways to assess the uncertainty of the forecasts in cases in which no model misfit is evidenced. This can be performed using simulated datasets generated with the parametric model, as was done here: in cases in which observation noise was large, forecasting accuracy was found to be poor (irrespective of the approach used, Figure 4C). This last result differs from the ones of Perretti et al. (2013c). The discrepancy between their and the present results is due to the fact that forecasts were here grounded on observed population sizes (including observation noise) rather than true population sizes as they did. This choice of using observed population sizes is justified by the fact that only this piece of information is available in real settings.

The level of observation noise was found to be a key determinant of forecasting accuracy (Figure 4), recalling the common sense that good forecasts require good data with minimized observation error. In real case studies, one cannot know a priori whether the level of observation noise is sufficiently low to ground a reliable forecast. The benefit of parametric approaches is that they enable to assess the data quality with respect to a forecast objective. Indeed, such an assessment can be performed by looking at the forecast accuracy in simulated datasets based on the parametric model with fitted-to-data observation noises (Figure 4 for a θ value equal to 1). Still, users of parametric approaches will gain at collecting independent data on observation error, since this helps increase the reliability and accuracy of parametric inferences [START_REF] Ives | Estimating community stability and ecological interactions from timeseries data[END_REF].

Proponents of the nonparametric approach argued that it provides better forecasts than parametric approaches in an array of noised and chaotic case studies (Perretti et al. 2013a,b).

They further argued that their approach was easier to apply than computationally demanding parametric approaches. The results of this study do not contradict these two arguments.

Nonparametric forecasts were indeed found to outperform parametric ones in a number of -but not all -cases with strongly reduced computational costs (Figure 4). However, no methods are yet available to assess the reliability of nonparametric forecasts which is likely to vary greatly (Figure 4). Not knowing when a forecast is likely to significantly fail is a major problem that should not be disregarded [START_REF] May | Complex systems: ecology for bankers[END_REF].

Circumventing the drawbacks of parametric forecasting

Parametric approaches do have real limitations. In particular, parametric inference can be challenging for complex models. A particular ad-hoc solution to this problem has been used in this study [START_REF] Pisarenko | Statistical methods of parameter estimation for deterministically chaotic time series[END_REF]. But many other ways to tackle parametric inference for complex models have been developed in recent years, notably methods based on computerintensive simulations like approximate Bayesian computation (Beaumont 2010, Jabot and[START_REF] Jabot | Analyzing tropical forest tree species abundance distributions using a nonneutral model and through approximate Bayesian inference[END_REF] or synthetic likelihoods [START_REF] Wood | Statistical inference for noisy nonlinear ecological dynamic systems[END_REF]. These methodological innovations are further eased with dedicated statistical packages (e.g., [START_REF] Csilléry | abc: an R package for approximate Bayesian computation (ABC)[END_REF][START_REF] Jabot | EasyABC: performing efficient approximate Bayesian computation sampling schemes using R[END_REF]) and should drastically reduce limitations to parametric inference. Besides inference, model development has been eased in recent years with progress in techniques of sensitivity analysis [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF]) and more generally by the synthesis of good modeling practices [START_REF] Grimm | Towards better modelling and decision support: documenting model development, testing, and analysis using TRACE[END_REF].

Valuing the advantages of nonparametric methods

The demonstration by Perretti et al. (2013a,b,c) of the great potential of nonparametric approaches suggests that a fruitful avenue for future research might be to combine the respective strengths of both parametric and nonparametric approaches. For instance, one could jointly use parametric and nonparametric approaches, using the first to assess forecast reliability, and the second to obtain more efficient forecasts in cases they are outcompeting parametric ones.

Nonparametric embedding methodologies may also be useful for model development, by helping choosing the appropriate time-delays to use in process-based approaches. Another line of research might be to include an observation process in the nonparametric S-map methodology.

This may enable to refine the recovery of the chaotic attractor in the embedded time series, and may enable to assess the forecast accuracy of the nonparametric approach. This last perspective should be seen as the key priority for proponents of the nonparametric approach, since the lack of methods for forecast reliability assessment is currently the Achilles' eel of the approach.

  computed, and divided by the standard deviation of the reference time series to obtain standardized root mean square errors (SRMSE, see Hartig and Dormann 2013 for details, and the Appendix for the R scripts used). Constant forecasts equal to the population temporal mean would lead to a SRMSE equal to one. Therefore, forecasts with a SRMSE larger than one can be considered to have failed. Note that, contrary to what Perretti et al. (2013b,c) and Hartig and Dormann (2013) did, predictive time series were not launched from the true population sizes X(t), but from a maximum-likelihood estimator of the current population size. This choice ensures a more realistic setting where only past observed population sizes N(t) are known. Also note that it would have been possible to use as starting values for the validation time series the sampled MCMC values for the true population sizes X(t), but it would have implied to repeat the parametric fitting for each starting value of the validation time series, which would have been here prohibitively costly computationally.

Figure 1 .

 1 Figure 1. Histograms of model checking p min values when the parametric fit is based on the "true" model in white (data and model fitted based on θ=1), or on a "wrong" model in grey (data generated with θ=4 while the model fitted uses θ=1). Panels A, B and C correspond respectively to cases of low, middle and large observation noise s obs .

Figure 2 .

 2 Figure 2. Quantile values of the model checking statistics p min as a function of the discrepancy between the model fitted (θ=1) and data (simulated with θ). The boxes and middle lines correspond to the 25, 50 and 75 percentiles of the p min distributions. Panels A, B and C correspond respectively to cases of low, middle and large observation noise s obs .

Figure 3 .

 3 Figure 3. Forecast errors measured by the standardized root mean squared error (SRMSE) as a function of the number of time steps between current data and the predictions, using either the parametric (solid lines) or the nonparametric approach (dotted lines). Vertical bars stand for the 95% confidence interval of the SRMSE. Panels A, B, C (D, E, F) correspond to the low (large) noise cases. Panels A and D (respectively B and E, C and F) correspond to data generated with θ equal to 1, 1.2 and 4.

Figure 4 .

 4 Figure 4. Forecasting efficiency measured by Area SRMSE (see methods) as a function of the discrepancy between the model fitted (θ=1) and data (simulated with θ). Panels A, B and C correspond respectively to cases of low, middle and large observation noise s obs . Crosses and plain lines stand for parametric forecasts. Circles and dashed lines stand for nonparametric forecasts.

  and 4). For real case studies, it is difficult to advise the use of a particular threshold below which a model should be considered to fail at fitting the data. Such procedures enable to detect model misfit, when a posterior predictive p-value is clearly lower than the distribution of posterior predictive p-values obtained on data generated with the parametric model (e.g., Figure2A and B, for θ values equal or larger than 1.5). In such cases, they will indicate particular patterns in the data that the model may not reproduce accurately, and therefore guide model improvements. Model checking procedures also confirm a good model fit to data, when the posterior predictive p-value is within the central range of the reference distribution. In less clearcut cases, model checking still provides ecologists with a quantification of the model fit quality.

  . Interpreting them as standard p-values would thus lead to too liberal model checking procedures. A second difficulty arises here since five p-values are combined in p min . Therefore, model checking statistics are not directly interpreted. Instead, each posterior predictive p-value is compared to the ones obtained on data sets generated with the model used for the fit. If the posterior predictive p-value obtained with the real data is lower than those obtained on data sets generated with the model, this will indicate a poor model fit quality of the real data.

Acknowledgements

FJ was funded by the Irstea INDECO project DynIndic. I thank Florian Hartig, the editor and an anonymous reviewer for their insightful comments on a previous version of this manuscript.


Thibaut Fréjaville
email: thibaut.frejaville@gmail.com

Thomas Curt

Christopher Carcaillet

Tree cover and seasonal precipitation drive understorey flammability in alpine mountain forests

Keywords: fire intensity, flammability traits, FlamMap model, mountain forest fuel, pyrogeography, redundancy analysis

Aim

Little is known about the understorey flammability of European mountain forests. We aimed to disentangle the relative effects of climate, vegetation structure and composition on the fuel-driven variation in fire spread and intensity.

Location The western Alps.

Methods Fire spread and intensity were simulated under constant moisture and weather conditions for a wide range of understorey fuel parameters measured in the litter, grass and shrub layers. Simulation outputs were used to compare understorey flammability between forest ecosystem types (FET). FETs were prior discriminated from a co-inertia analysis between composition and the environment (vegetation structure and climate). Relationships between these factors, fuel properties and understorey flammability were then tested using partial regression analyses.

Results

The most flammable forests display an open canopy ('dry-subalpine' and 'open-Mediterranean') and grow in areas with dry autumns, wet and cold springs. Fire spread and intensity are controlled by the trade-off between tree cover and the load of dead (litter) and live biomass (grass and shrub). Otherwise, fire intensity increases with precipitation seasonality (magnitude between seasons); rainy springs enhance biomass growth whereas dry climates (especially in autumn) promote shrub biomass and stimulate litter accumulation and residence (higher litterfall and lower decomposition). Interestingly, we found a positive relationship between fire intensity and the proportion of conifers in the community that disappeared after accounting for tree cover, highlighting that, in the Alps, the open canopy structure of needleleaved forests makes them more flammable than the broadleaved due to a higher amount and continuity of surface fuels.

Main conclusions

Interrelationships between tree cover, precipitation seasonality and species composition govern the understorey flammability of mountain forests. We add further evidences that tree cover strongly constrains fire spread by driving the amount and the type of surface fuels, what suggests strong influences of land-use changes on flammability patterns.

INTRODUCTION

Plant flammability is an important fire driver with a fundamental role in terrestrial ecosystems [START_REF] Bond | Fire as a global "herbivore": the ecology and evolution of flammable ecosystems[END_REF]. It , and is related to morphological traits of plants, including canopy architecture, fine fuel biomass and surface area-to-volume ratio (Bowman et al., 2014a). Flammability traits differ between plants [START_REF] Papió | Comparative study of the aerial structure of five shrubs of Mediterranean shrublands[END_REF][START_REF] Behm | Flammability of native understorey species in pine flatwood and hardwood hammock ecosystems and implications for the wildland-urban interface[END_REF][START_REF] Fréjaville | Bark flammability as a fire-response trait for subalpine trees[END_REF] and affect the spread and the intensity of fires [START_REF] Santana | Fuel structural traits modulating soil temperatures in different species patches of Mediterranean Basin shrublands[END_REF][START_REF] Schwilk | Scaling from leaf traits to fire behaviour: community composition predicts fire severity in a temperate forest[END_REF]. Consequently, certain plant communities are more likely to burn than others [START_REF] Bond | Fire and Plants[END_REF]. The composition of plant communities substantially modifies fire regime [START_REF] Engber | Patterns of flammability of the California oaks: the role of leaf traits[END_REF] and, reciprocally, fire regime alters the community flammability [START_REF] Curt | Litter flammability in oak woodlands and shrublands of southeastern France[END_REF][START_REF] Ganteaume | Effects of vegetation type and fire regime on flammability of undisturbed litter in Southeastern France[END_REF]. The role of strong contrasts of vegetation flammability in determining past and future fire activity has recently been emphasized [START_REF] Girardin | Vegetation limits the impact of a warm climate on boreal wildfires[END_REF][START_REF] Paritsis | Habitat distribution modeling reveals vegetation flammability and land use as drivers of wildfire in SW Patagonia[END_REF]. However, predicting spatial or temporal changes in vegetation flammability still needs to understand how flammability traits vary with the environmental conditions (e.g. climate). Identifying the environmental drivers of plant flammability traits, or more generally modelling the biogeographical variation of plant functional traits, appears thus critical to predict ecosystem responses to ongoing global changes [START_REF] Lavorel | Predicting changes in community composition and ecosystem functioning from plant traits: revisiting the Holy Grail[END_REF][START_REF] Reichstein | Linking plant and ecosystem functional biogeography[END_REF].

Because wildfire propagates through multispecific vegetation, the plant community appears as the most appropriate unit for flammability analysis (Bowman et al., 2014a). At the community level, flammability can be assessed by physical characteristics of fire behavior, i.e., the rate of fire spread and the intensity [START_REF] Anderson | Forest fuel ignitibility[END_REF][START_REF] Gill | Flammability of Australian forests[END_REF], which enable scaling up, comparison and linkage with fire behavior models [START_REF] Fernandes | Plant flammability experiments offer limited insight into vegetation-fire dynamics interactions[END_REF]. Indeed, community-level differences in flammability reflect the combined effects of varying sets of flammability traits carried (enhanced?) by living and dead plant material of the constituent species. The flammability of a community is therefore controlled by intrinsic fuel properties of plant assemblages, which determine the spread and intensity of fire under given conditions of topography and weather. While short-term variations in temperatures, humidity or wind govern the fuel moisture and ignitability [START_REF] Flannigan | Implications of changing climate for global wildland fire[END_REF], mid-to long-term climate should constrain fire spread and intensity by controlling biomass growth and fuel spatial arrangement (vegetation structure, [START_REF] Pausas | Fuel shapes the fire-climate relationship: evidence from Mediterranean ecosystems[END_REF] and by constraining the composition of plant communities [START_REF] Bradstock | A biogeographic model of fire regimes in Australia: current and future implications[END_REF]. Consequently, disentangling how climate acts on the overall community flammability is a challenging task [START_REF] Parisien | Environmental controls on the distribution of wildfire at multiple spatial scales[END_REF].

In mountains, a strong increase in fire activity is expected in response to climate warming [START_REF] Westerling | Warming and earlier spring increase western U.S. forest wildfire activity[END_REF][START_REF] Moritz | Climate change and disruptions to global fire activity[END_REF], notably in the western Alps where fires become more frequent and large [START_REF] Fréjaville | Spatiotemporal patterns of changes in fire regime and climate: defining the pyroclimates of south-eastern France (Mediterranean Basin)[END_REF]. In these ecosystems, understanding the environmental drivers of forest flammability along the wide climate gradients and the diverse vegetation communities will improve our ability to manage fire in mountains. To address it, we aimed i) to compare understorey fuel properties between forest ecosystems in the Alps and, ii) to disentangle the relative influence of community composition and the environment (forest structure and climate factors) on understorey flammability and on the underlying fuel structural traits, and iii) to determine the main flammability drivers. We applied the following framework (Fig. 1) in the western Alps. First, the main variability in community composition, forest structure and climate was identified by defining forest ecosystem types (FET). Second, based on Rothermel's equations (1972), surface fuel parameters were collected under a stratified sampling scheme within FETs to simulate fire behavior along the range of habitat conditions and elevation. Then, understorey flammability was calculated by the trend of surface fuels to burn fast and with high intensity under given external conditions of temperature, humidity, wind, aspect and slope. Finally, differences in understorey flammability between forest ecosystems were related to community composition and to the environmental conditions of forest structure and climate using partial regression analysis. As flammability was calculated for a given weather, testing the influence of climate parameters on flammability relies here on mid-to long-term variation in fuel characteristics, i.e., composition, load and arrangement) other than moisture (short-term variation). We focused on understorey flammability, because fire regimes were typified by low-intensity surface fires in the Alps [START_REF] Genries | The function of surface fires in the dynamics and structure of a formerly grazed old subalpine forest[END_REF] as in other mountains [START_REF] Schoennagel | The interaction of fire, fuels, and climate across Rocky Mountain forests[END_REF][START_REF] Rother | Climatic influences on fire regimes in ponderosa pine forests of the Zuni Mountains, NM, USA[END_REF], highlighting the dominant function of surface fuels (litter, herbs and shrubs) for fire spread. We hypothesize that (H1) the wide array of climates and vegetation in the Alps results in contrasting fuels properties between forest ecosystems (in terms of quality and spatial arrangement), which leads to high variation in the spread and intensity of fires. Further we hypothesize that (H2) this flammability variation could be predicted from vegetation and climate factors that act on the fuel build-up or the community composition.

METHODS

Data sources

Structure and composition parameters of mountain forests in the western Alps were obtained from seven annual campaigns (2005)(2006)(2007)(2008)(2009)(2010)(2011) of the French national forest inventory (Inventaire forestier national, hereafter IFN). The IFN comprises a network of temporary plots established on a grid of c. 500×500 m. We focus on a study area of 31 710 km 2 (43°30'-46°24'N, 4°60'-7°40'E) that extends from the northern to the southern range of the French Alps (Fig. 2), within which the climate of the lowlands varies from warm-dry (Mediterranean) to cool-wet (oceanic) through continental. We excluded plots if any evidence of recent (<5 years) natural or anthropogenic disturbances such as logging, fire or windthrow. Further, we removed plots below 700 m a.s.l. and those from valley floor to focus on forests Author-produced version of the article published in Journal of Biogeography, 2016, 1, 1-12. The original publication is available at http://onlinelibrary.wiley.com/doi/10.1111/jbi.12745/full DOI: 10.1111/jbi.12745 situated on slopes. Parameters of climate, forest structure and abundance of woody species were analyzed from 2609 IFN plots evenly distributed across the French Alps (Fig. 2). IFN data include quantitative estimates of plant species abundances by phytosociological scores (Braun-Blanquet, 1961), horizontal covers of vegetation layers (herb, shrub, tree), and morphological variables such as tree diameter (at breast height) and height.

The forest composition (FC) is described here by the relative abundance of shrub and tree species (species composition) and of plant functional types (functional composition), i.e., based on needleleaf vs broadleaf, deciduous vs evergreen, trees vs shrubs. Prior to analyses, abundance data was transformed from Braun-Blanquet scores recorded in the IFN plots to median values of cover classes. By this way we ensured compatibility with Euclideandistance-based ordination methods [START_REF] Podani | Braun-Blanquet's legacy and data analysis in vegetation science[END_REF]. The forest structure component (FS) is described by the canopy characteristics (tree cover, tree density, basal area and mean tree height) and by the horizontal cover of shrubs and herbs. Climatic (CL) parameters were obtained from data downscaled at a 1×1 km resolution using a regional regression procedure and a 500×500 m digital elevation model [START_REF] Fréjaville | Spatiotemporal patterns of changes in fire regime and climate: defining the pyroclimates of south-eastern France (Mediterranean Basin)[END_REF]. Specifically, we used long-term averages of monthly precipitation, monthly temperatures and the annual degree-day sum over the growing season (values above 5°C, GDD5). We also used long-term averages (annual mean, monthly mean and annual 95th percentile) of the drought indices of the Canadian fire weather index system (van [START_REF] Wagner | Development and structure of the Canadian Forest Fire Weather Index System[END_REF]: the fine fuel moisture code (FFMC), the duff moisture code (DMC) and the drought code (DC). Climate data computed over the latest five years (2005)(2006)(2007)(2008)(2009) was preferentially used in analyzing the flammability drivers, which better captures climate-fuel relationships than longer series.

Ecosystems study

To insure that our flammability analysis was extended to the variability in climate and vegetation found in Alpine forests, understorey fuel parameters were sampled after different habitats were defined. Precisely, a broad classification of mountain forest ecosystems (FET) was carried out using a co-inertia analysis (OMI, [START_REF] Dolédec | Niche separation in community analysis: a new method[END_REF] of woody species assemblages (shrubs and trees) and a principal component analysis of the environment, i.e. climate and forest structure (see Appendix S1 in Supporting Information). It resulted in four FETs. The 'dry-subalpine' ecosystems (FET-1), which cover the inner continental zone of the Alps (Fig. 2), are dominated by Larix decidua, Pinus uncinata or Pinus cembra forming a tall and open canopy over a high cover of grass or dwarf shrubs, mainly Juniperus sibirica, Vaccinium spp. and Rhododendron ferrugineum. The 'moist-montane' ecosystems (FET-2), which cover the wet northern range (Fig. 2), are dominated by Abies alba, Picea abies or Author-produced version of the article published in Journal of Biogeography, 2016, 1, 1-12. The original publication is available at http://onlinelibrary.wiley.com/doi/10.1111/jbi.12745/full DOI: 10.1111/jbi.12745

Fagus sylvatica forming tall closed forests with a scarce understorey. At lower elevations, two other types gather low-stature forests of Quercus pubescens, Pinus nigra or Pinus sylvestris that uncover the warm and dry southern range under Mediterranean influences. The 'open-Mediterranean' ecosystems (FET-3) display open forests with moderate grass cover and high cover of small shrubs such as Juniperus communis, Genista cinerea and Amelanchier ovalis. The 'closed-Mediterranean' ecosystems (FET-4) have a higher tree cover with a more or less dense understorey of tall shrubs such as Buxus sempervirens and Acer species.

Fuel structural traits

Surface fuel parameters were then measured in 96 'Fuel' plots (Fig. 2) from a stratified random sampling design within FET, using as criterions the proximity to IFN plots of the FETs and favoring multiple locations along elevation gradients. Fuel measurement was based on standard procedures (e.g. [START_REF] Bessie | The relative importance of fuels and weather on fire behaviour in subalpine forests[END_REF] by following methods of [START_REF] Mcrae | Measurement and description of fuels and fire behaviour on prescribed burns: A handbook[END_REF] and [START_REF] Brown | Handbook for inventorying surface fuels and biomass in the Interior West[END_REF] along an equilateral triangle of 20 m aside. Downed woody debris were measured by the line intersect technique ( [START_REF] Wagner | The line intersect method in forest fuel sampling[END_REF]). Cover, height and load of herb and litter layers (i.e. the superficial fraction of fuel particles ø< 6 mm not yet decayed) were measured in nine quadrats of 1-m 2 evenly spaced along the sides of the triangle. Surface to volume ratio (SVr) was measured for litter and shrub leaves on the basis of geometric dimensions of leaf middle cross-area [START_REF] Hachmi | Alternative equations to estimate the surface-to-volume ratio of different forest fuel particles[END_REF]. SVr measurements were averaged from five leaves randomly sampled by forest plot for each dominant species (absolute cover >10%). Species composition and heights of shrubs (crown-base and -top heights) were measured every 2-m along the 60-m transect (triangle sides) for all individuals (fuel ladders) that vertically intercepted the transect.

In the case of a multi-layered shrub stratum, considering all shrubs would distort the 'realism' of fire behavior simulations because if fire is unable to bridge gaps between strata, then the higher fuels are unavailable. Therefore, we developed a procedure to detect the vertical arrangement of shrubs, i.e. to detect whether shrub individual heights follow a multi-modal distribution using kernel density functions (Appendix S2). Fire simulations were conducted in each detected vertical arrangement of shrub layers. The belonging of each shrub layer to the fuel bed was assessed from the critical fire intensity criterion of [START_REF] Wagner | Conditions for the start and spread of crown fire[END_REF] and iterative fire behaviour simulations (Appendix S2). Shrub load was prior estimated from the fine fuel fraction (leaves and twigs ø < 6 mm) and the crown volume of dominant species:

L Sx =0 .25⋅ ∑ j=1 S BD j⋅¿ phyV j ¿ (1)
Author-produced version of the article published in Journal of Biogeography, 2016, 1, 1-12. The original publication is available at http://onlinelibrary.wiley.com/doi/10.1111/jbi.12745/full DOI: 10.1111/jbi.12745

where LSx is the load of the shrub (multi-)layer x, BDj the bulk density and phyVj the cubicshape phyto-volume of the constituent dominant species j, and 0.25 a correction coefficient (Appendix S2). Community-level values of SVr were computed by weighting the mean value of the constituent species j by its relative abundance Rj:

SVr = ∑ j=1 S SVr j⋅¿ R j / ∑ j=1 S R j ¿ (2)
In the case of shrub leaves, Rj is the relative phyto-volume of the j th shrub species in the (multi-)layer, whereas in the case of tree litter leaves, Rj is the relative basal area of the j th tree species in the community.

Community flammability quantification

Anderson's definition (1970) of flammability encompasses three components describing the fuel ability to ignite (ignitability), the intensity of combustion (combustibility) and its sustainability. Fire line intensity (FLI in kilowatt per meter, kW.m -1 ) reflects combustibility while the concept of ignitability is obvious in fire rate of spread, ROS, expressed in meter per minute, m.min -1 [START_REF] Gill | Flammability of Australian forests[END_REF], which idealized successive ignitions of fuel fragments from surrounding burning ones. The heat released per unit area, HPA, in kilojoule per square meter, kJ.m -2 , expresses both combustibility and sustainability because the amount of heat released depends on both intensity and duration of combustion.

Understorey flammability was quantified by the multivariate variation of ROS, FLI and HPA. Fire behavior was simulated on one virtual terrain (slope = 30%, wind blowing uphill) for different weather scenarios of dead fuel moisture and wind speed using the FlamMap program [START_REF] Finney | An overview of FlamMap fire modeling capabilities[END_REF]. Canopy characteristics (tree cover, mean height) were used to take into account the buffering effect on surface wind speed (i.e. the wind adjustment factor; [START_REF] Andrews | Modeling wind adjustment factor and midflame wind speed for Rothermel's surface fire spread model[END_REF]. Live fuel moisture was set at 100%. Changes in dead fuel moisture and wind speed did not significantly change neither the relative differences in fire behavior between FET, nor the relationships with environmental conditions (data not shown). Consequently, we focus here on an intermediate severity scenario for the sake of conciseness, i.e., under a 30 km.h -1 wind speed and 7% dead fuel moisture content. Fire behavior simulations were performed independently for all fuel models, i.e. for one to several homogeneous fuel complexes per community (Appendix S2). When several fuel models were designed in a community, we averaged fire simulation outputs according to the relative cover of the fuel model in the community.

Author-produced version of the article published in Journal of Biogeography, 2016, 1, 1-12. The original publication is available at http://onlinelibrary.wiley.com/doi/10.1111/jbi.12745/full DOI: 10.1111/jbi.12745 montane, FET-2) depict lower amount of live fine fuels, i.e. herb and shrub within the fuel bed (p < 0.05). Dwarf shrubs of dry-subalpine ecosystems have a higher bulk density and tree litter leaves a higher SVr than other FETs (p < 0.05).

Considering the combined effect of these fuel structural traits, fire simulations indicated an elevated flammability for the open forests (FET-1 dry-subalpine and FET-3 open-Mediterranean) that exhibited the highest fire rate of spread (ROS), fireline intensity (FLI) and heat released per unit of area (HPA) for given fuel moistures and wind speed (Fig. 3). Closed-Mediterranean ecosystems (FET-4) are characterized by an intermediate flammability with a lower mean ROS (Fig. 3a), but similar FLI and HPA values than open forest ecosystems (Fig. 3bc). The poor-understorey moist-montane ecosystems (FET-2) are the less flammable.

Environmental drivers of understorey flammability

Boosted redundancy analysis indicated that parameters of forest composition, structure or climate explain most of the flammability variance (quantified by the multivariate variation in ROS, FLI and HPA) within and between FETs (R²mean = 0.62, pmean < 0.001, Fig. 4a). The interaction between these three components explained the main part of flammability variance (R²mean = 0.39) compared to the pure fractions (variance that is exclusively explained by one component). Nonetheless, differences in flammability are also exclusively related to forest structure (R²mean = 0.11, pmean < 0.05, Fig. 4a), and to a lesser extent to forest composition or climate (pure fractions: R²mean = 0.07 and 0.04 respectively, pmean ≤ 0.05, Fig. 4a). The interactions between forest structure, composition and climate explain the main part of fuel parameters variance (Table 2), except for the fuel load (litter, herb and shrub) that varied mainly with climate (half of explained variance, p < 0.01, Table 2). Community composition drove most of the variance in surface area to volume ratio of shrub leaves (71% of explained variance, p < 0.05, Table 2).

Stepwise regressions revealed that tree cover, precipitation seasonality (annual to spring precipitation ratio, or the magnitude between seasons) and forest functional composition chiefly explain the variance of flammability between forests (adjusted R² = 0.56, p < 0.001, Fig. 4b). It is better explained by the interaction between these three factors (23%) and by canopy opening (variance that is exclusively explained by tree cover: adjusted R² = 0.22, p < 0.001, Fig. 4b). In particular, FLI decreases non-linearly with the tree cover by increasing until ca. 40% of tree cover threshold and then decreasing with canopy closure (r² = 0.44, p < 0.001, Fig. 5a). A multi-linear regression indicated that FLI increases also with the precipitation seasonality (t = 4.5, p < 0.001, Fig. 5b

), showing that open canopy forests are

Author-produced version of the article published in Journal of Biogeography, 2016, 1, 1-12. The original publication is available at http://onlinelibrary.wiley.com/doi/10.1111/jbi.12745/full DOI: 10.1111/jbi.12745 of ca. 40%, above which fire intensity decreased (Fig. 5a). This threshold was already detected in relation to burned areas in southern Africa [START_REF] Archibald | What limits fire? An examination of drivers of burnt area in Southern Africa[END_REF]. This was supported by a nonlinear response of grass productivity to tree cover that inflects around 35-40% cover [START_REF] Scholes | Convex relationships in ecosystems containing mixtures of trees and grass[END_REF] in reaction to decreasing light intensity and increasing rainfall interception by deep litter layers [START_REF] Anderson | Herbaceous response to canopy cover, light intensity, and throughfall precipitation in coniferous forests[END_REF]. In the western Alps, a trade-off between the amount of litter and surface fuels (herb and shrub) might optimize the fire intensity at an intermediate canopy cover (Fig. S2.2), i.e.,, for a high amount of grass and shrub fuel on a continuous litter bed. This supports observations in savannas, where tree cover decreased the probability of burning and fire intensity by limiting the availability of fine fuels [START_REF] Trauernicht | Tree cover-fire interactions promote the persistence of a fire-sensitive conifer in a highly flammable savanna[END_REF], in addition to the sheltering effect of canopy closure on microclimatic conditions [START_REF] Ray | Micrometeorological and canopy controls of fire susceptibility in a forested Amazon landscape[END_REF][START_REF] Hoffmann | Fuels or microclimate? Understanding the drivers of fire feedbacks at savanna-forest boundaries[END_REF].

To test the relationships between tree cover and understorey flammability that we found using fire simulations, we analysed the distribution of lightning fire sizes for the period 1980-2011 according to the tree cover of surrounding unburned IFN plots (Appendix S3). Among Mediterranean mountain forests of the Southern Alps, we found a higher burned area ratio in open (FET-3) than in closed (FET 4) ecosystems (Appendix S3), comforting the simulated differences of fire behaviour between them (Fig. 3). Further, we found that forests display a more open canopy (median tree cover value of 40%) near the largest lightning fires (> 10 ha) compared to areas where lightning fires were small or absent (median tree cover value of 80%; X²df =2 = 9.37, p = 0.009, Appendix S3). These independent validation tests strengthen the relationships we found above between tree cover, understorey flammability and surface fire spread in the Alps (Fig. 5, Fig. S2.2).

Further, fire intensity varied between forests mainly due to interactions between tree cover, precipitation seasonality and functional composition (Fig. 4). For instance, we demonstrated that coniferous (needleleaved) forests are more flammable than broadleaved deciduous forests, a widely accepted assumption [START_REF] Bond | Fire and Plants[END_REF][START_REF] Girardin | Vegetation limits the impact of a warm climate on boreal wildfires[END_REF], which is however significantly related to differences in tree cover (Fig. 7). This observation highlights the ability of broadleaved deciduous trees to form more closed canopies and therefore less flammable forests. We also found an increase in the amount of fine fuels (litter, herb and shrub) with precipitation seasonality, i.e. with increasing spring rainfall and decreasing autumn rainfall, leading to more intense fires in climates characterized by wet springs and dry autumns, where open forest ecosystems dominate (FET-1 and FET-3, higher position on precipitation seasonality gradient, Fig. 6). This supports the principle that climate seasonality controls the flammability of mountain forests by promoting open canopies, biomass growth under rainy springs and necromass accumulation under dry autumns. Indeed, we found a negative correlation between the amount of fine litter, autumn rainfall and spring temperatures, suggesting a decrease in litter decomposition rate towards drier or colder areas [START_REF] Meentemeyer | World patterns and amounts of terrestrial plant litter production[END_REF][START_REF] Aerts | Climate, leaf litter chemistry and leaf litter decomposition in terrestrial ecosystems: a triangular relationship[END_REF][START_REF] Keane | Biophysical controls on surface fuel litterfall and decomposition in the northern Rocky Mountains, USA[END_REF]. In a lesser extent, litter load varied Author-produced version of the article published in Journal of Biogeography, 2016, 1, 1-12. The original publication is available at http://onlinelibrary.wiley.com/doi/10.1111/jbi.12745/full DOI: 10.1111/jbi.12745 according to tree species composition and functional traits (p < 0.05, Table 2), as already found in the northern Rocky Mountains [START_REF] Keane | Biophysical controls on surface fuel litterfall and decomposition in the northern Rocky Mountains, USA[END_REF] and in dry-subalpine ecosystems of the Alps [START_REF] Blarquez | Needle accumulation rate modelbased reconstruction of palaeo-tree biomass in the western subalpine Alps[END_REF]. Therefore we add strong evidences that community flammability in surface-fire ecosystems is controlled by the canopy cover and the interplay with climate (precipitation seasonality), fuel connectivity and litter bed characteristics (Bowman et al., 2014a).

Otherwise, we found that fire intensity increased with drought (independently to fuel moisture), especially in autumn (Table S2.2), but following an unimodal relationship with the most extreme annual conditions (95 th percentile of DC, Fig. 6c). These findings highlight that dry conditions enhance understorey flammability in the Alps by promoting shrub biomass instead of tees biomass, and litter accumulation (Table S2.1) whereas fuel limitation by drought may also constrain fire in supra-Mediterranean ecosystems. Under a given fuel moisture, we also found that forest ecosystems inhabiting the colder climates (low annual, winter and spring temperatures) burn more readily. These findings illustrate that despite temperature is recognized as the most important variable affecting wildfire, with warmer temperatures leading to increased fire activity [START_REF] Flannigan | Implications of changing climate for global wildland fire[END_REF], a warmer climate could also inhibit fire by promoting less flammable broadleaf [START_REF] Girardin | Vegetation limits the impact of a warm climate on boreal wildfires[END_REF]. For instance, surface-area-to-volume ratio of litter leaves decreased with winter and spring temperatures (Table S2.1), because composition-driven changes in this flammability trait (Table 2) mostly relate to these climatic gradients.

By controlling fuel dryness, biomass growth, decomposition of dead plant material and the community composition, climate change should strongly alter fire spread, intensity and ultimately fire regime in mountain ecosystems by promoting synergistic trends (e.g., higher fuel load with drier conditions) versus antagonist ones (e.g., temperature-or drought-induced fuel limitation). Otherwise, the open canopy structure of the most flammable forests likely mirrors transient afforestation processes of abandoned pastures and eventually ploughed areas at low (open-Mediterranean ecosystems) and high (dry-subalpine ecosystems) elevations [START_REF] Tasser | Impact of land use changes on mountain vegetation[END_REF][START_REF] Chauchard | Patterns of land-use abandonment control tree-recruitment and forest dynamics in Mediterranean mountains[END_REF]. Therefore, past land-use changes in European mountain forests might have governed current flammability patterns at the landscape scale [START_REF] Zumbrunnen | Human impacts on fire occurrence: a case study of hundred years of forest fires in a dry alpine valley in Switzerland[END_REF]. Consequently, the human-induced ongoing global changes will likely impact the occurrence, spread and intensity of fires in mountain forests by promoting drier conditions [START_REF] Zumbrunnen | Linking forest fire regimes and climate -A historical analysis in a dry inner Alpine valley[END_REF][START_REF] Fréjaville | Spatiotemporal patterns of changes in fire regime and climate: defining the pyroclimates of south-eastern France (Mediterranean Basin)[END_REF] and by inducing changes in the main environmental drivers of understorey flammability, i.e. the tree cover, the precipitation seasonality, and the relative abundance of broadleaved versus needleleaved trees. These findings also imply that under similar weather and fuel moisture conditions (as simulated here), fire spread and intensity should greatly differ across Alpine forests, in a predictable manner. We found a higher fire susceptibility in open forests of southern Mediterranean mountains and, more surprisingly, in the subalpine ecosystems of the Tables Table 1: Fuel parameters (mean ±sd) of mountain Forest Ecosystem Types (FET) of the western Alps. Between-FET comparisons were performed from ANOVA and HSD post-hoc tests; significant (p < 0.05) F and p-values are highlighted in bold with differences reported by different letters. Fuel parameters were fourth-root transformed before comparison tests. Fuel parameters involved here the set of dead and live fine biomass available for burning by the surface fire front (i.e. within the fuel bed, see the Methods). Mean annual temperatures and precipitation of FET were also presented. SVr, surface area to volume ratio

Figure captions Figure 1 :

 captions1 Figure captionsFigure1: Analytical framework representing the approach used to test differences in understorey flammability between forest ecosystems according to plant composition and environmental factors (forest structure and climate).

Figure 2 :

 2 Figure 2: Map of Forest Ecosystem Types (FET) in the western Alps: location of IFN (+) and fuel sampling (○) sites (left), and contour of FET (right).

Figure 3 :

 3 Figure 3: Flammability of Forest Ecosystem Types (FET). Weather-constant simulation outputs (n = 96) of fire rate of spread (ROS, a), fire line intensity (FLI, b) and heat released per unit area (HPA, c) were compared between FET by ANOVA and HSD post-hoc tests. Differences (p < 0.05) are reported by different letters. Note a log-scale for vertical axes.

Figure 4 :

 4 Figure4: Variance partitioning of community flammability (matrix of log-transformed ROS, FLI and HPA) into unique and common ('inter.') contributions of (a) the sets of predictors (forest structure 'FS', composition 'FC', and climate 'CL') or (b) only the most parsimonious ones. Variance was partitioned using partial redundancy analyses, with a bootstrap procedure in (a) from which bootstrapped values (n = 50) of adjusted R² were compared between unique and common contributions by ANOVA and HSD post-hoc tests, with significant differences (p < 0.05) reported by different letters.

Figure 5 :

 5 Figure5: Relationships between fire-line intensity (log-transformed) and tree cover (a) with precipitation seasonality (spring to annual precipitation ratio) as a covariable (b). Variation of fire-line intensity with tree cover was regressed by a two-order polynomial function. Means (±sd) of FET are depicted: 'dry-subalpine' (FET-1) in black, 'moist-montane' (FET-2) in red, 'open-Mediterranean' (FET-3) in green, 'closed-Mediterranean' (FET-3) in blue.

Figure 6 :Figure 7 :

 67 Figure6: Partial relationships between fire line intensity (FLI, log-transformed) and precipitation seasonality (spring to annual precipitation ratio, a), autumn precipitation (b) and annual 95th percentile of daily DC (c). Partial relationships between these climate gradients and the fine fuel load were also shown (bottom). The part of variance explained by the environmental components X1 and X2 independent to the studied variable x (X1+X2 | x) was removed from redundancy analyses before modelling partial relationship on the unexplained variance fraction (indicated in brackets). Means (±sd) of FET are depicted: 'dry-subalpine' (FET-1) in black, 'moist-montane' (FET-2) in red, 'open-Mediterranean' (FET-3) in green,

Figure 1 Figure 3 Figure 4

 134 Figure 1

Table 2 :

 2 Author-produced version of the article published in Journal of Biogeography, 2016, 1, 1-12. The original publication is available at http://onlinelibrary.wiley.com/doi/10.1111/jbi.12745/full DOI: 10.1111/jbi.12745 Variance of forest fuel parameters (fourth-root-transformed) explained by forest structure, composition and climate components. These three environmental components were described by the ≤ 12 most parsimonious variables for each. Total explained variance was decomposed from 50 (partial) redundancy analyses into pure and interaction fractions. The pure fraction describes the amount of variance that is exclusively explained by one component (i.e. after accounting for other components). Significance was based on Monte Carlo permutation tests. Variance and p-values were averaged. * p < 0.05, ** p < 0.01, *** p < 0.001. SVr, surface area to volume ratio

				FET-1	FET-2	FET-3	FET-4
				'dry-subalpine' 'moist-montane'	'open-Med'	'closed-Med'
		Unit	F3,92 (p-value)	n = 27	n = 16	n = 29	n = 24
	annual temperature	(°C)	70.3 (< 0.001)	4.5 (1.5) d	6.4 (1.2) c	8.3 (1.3) b	8.6 (1.2) a
	annual precipitation	(mm) 27.6 (< 0.001)	1193 (267) b	1551 (251) a	994 (112) d	1058 (133) c
	fine litter load		2.1 (0.104)	6.68 (3.49)	5 (1.91)	5.42 (2.84)	4.68 (2.18)
	dead woody load		1.8 (0.153)	10.46 (4.05)	11.42 (4.75)	8.75 (4.03)	9.36 (2.9)
	herb load shrub fine fuel load	(t.ha -1 )	8.2 (< 0.001) 4.4 (0.006)	0.52 (0.44) a 0.56 (0.69) a	0.35 (0.89) bc 0.02 (0.07) b	0.38 (0.43) ab 0.79 (1.18) a	0.16 (0.25) c 0.74 (1.16) ab
	live fine fuel load		6.6 (< 0.001)	1.08 (0.84) a	0.37 (0.9) b	1.4 (1.61) a	0.9 (1.2) ab
	total fuel load		1.0 (0.382)	11.54 (4.4)	11.79 (4.83)	10.15 (4.36)	10.26 (3.43)
	fuel depth	(cm)	3.3 (0.024)	29.66 (12.01) a	17.54 (12.86) b	36.59 (28.71) a	29.52 (26.42) ab
	shrub bulk density	(kg.m -3 ) 7.1 (< 0.001)	2.04 (0.96) a	0.98 (0.56) ab	0.88 (0.61) b	1.06 (0.72) b
	litter leave SVr shrub leave SVr	(cm -1 )	7.2 (< 0.001) 77.11 (28.16) a 2.8 (0.051) 65.77 (7.01)	67.69 (12.72) ab 76.32 (1.75)	54.31 (14.44) b 58.58 (10.85)	57.42 (18.63) b 64.17 (18.03)

ACKNOWLEDGEMENTS

Financial support was provided by the FUME Project under the European Union's Seventh Framework Programme (FP7/2007(FP7/ -2013) ) and by grants from the National Research Institute of Science and Technology for Environment and Agriculture. We acknowledge all persons who have contributed to the French national forest inventory database. We also thank Jean-Philippe Vidal for providing climate data and Jean-Michel Lopez, Roland Estève, Fabien Guerra, Christian Ripert and Annabelle Rivoal for their invaluable assistance in fieldwork.

Author-produced version of the article published in Journal of Biogeography, 2016, 1, 1-12. The original publication is available at http://onlinelibrary.wiley.com/doi/10.1111/jbi.12745/full DOI: 10.1111/jbi.12745

Variance partitioning

We first performed ANOVA and Tukey's HSD post-hoc tests to compare understorey flammability and fuel parameters between FET. To achieve normality and variance homogeneity, we applied a log-transformation of fire behavior metrics (ROS, FLI and HPA) and a fourth-root transformation of fuel parameters before comparison tests.

Then, the variance in understorey flammability was partitioned into unique and common contributions of the sets of predictors (climate, composition and forest structure) using partial redundancy analysis [START_REF] Peres-Neto | Variation partitioning of species data matrices: estimation and comparison of fractions[END_REF]. The flammability matrix included ROS, FLI and HPA community-level values (multivariate response). For each set of predictor (FC, FS, CL), we applied a forward selection of variables [START_REF] Blanchet | Forward selection of explanatory variables[END_REF]. Variance and pvalues were averaged from 50 partial redundancy analyses with Monte Carlo permutation tests, the residuals under the 'full' model being permuted 10 000 times [START_REF] Legendre | Numerical Ecology[END_REF], carried out on 70% of randomly selected sites (bootstrap procedure, see [START_REF] Cottenie | Integrating environmental and spatial processes in ecological community dynamics[END_REF]. Fuel parameters were modelled in the same way, expect that the procedure was performed separately for each (univariate response) and that the sets of predictors was limited to the most parsimonious subset of ≤12 variables each, using stepwise model selection (both backward and forward) with the Akaike Information Criterion (AIC).

From these redundancy analyses we also performed partial regressions to study to what extent and how one environmental factor drives changes in flammability and underlying fuel properties, i.e. by accounting for the influence of other drivers. Specifically, the effect of one predictor on one response variable (FLI or fuel load) was modelled by excluding the unique contributions of other confounding factors [START_REF] Peres-Neto | Variation partitioning of species data matrices: estimation and comparison of fractions[END_REF].

All analyses and calculations were realized with the R version 3.1.2 (R Core Team, 2013), using the ade4 and vegan packages for FET classification and variance partitioning.

RESULTS

Understorey flammability

Comparison of surface fuel parameters indicated significant biogeographic contrasts amongst FETs (Table 1). Litter is the most important contributor to the fuel bed by accounting for the higher proportion to total fine fuel load (Fig. S2.1). However, no significant differences in fine fuel load was detected between FET for litter, dead woody debris and thus total fuel (live and dead; p > 0.05, Table 1). In contrast, open FET (FET-1 dry-subalpine and FET-3 open-Mediterranean) show higher loads of herb (p < 0.05), whereas dense tall forests (moist-Author-produced version of the article published in Journal of Biogeography, 2016, 1, 1-12. The original publication is available at http://onlinelibrary.wiley.com/doi/10.1111/jbi.12745/full DOI: 10.1111/jbi.12745 the most flammable especially in areas where precipitation are abundant in spring and rare in other seasons (adjusted R² = 0.54, p < 0.001, Fig. 5b).

Among climatic parameters, partial regressions revealed that FLI increases first with precipitation seasonality (p < 0.001, Fig. 6a), then decreases with autumn precipitation (p < 0.001, Fig. 6b) and follows a humped relationship with the 95th percentile of the DC (p < 0.001, Fig. 6c). Total fine fuel load exhibits similar relationships with these climatic gradients (p < 0.001, Fig. 6def). Specifically, higher load of fine litter were found in areas with drier autumns and colder springs (p < 0.01, Table S2.1). Mean winter, spring and annual temperatures exhibit negative correlations with litter load, herb load, shrub bulk density and surface to volume ratio of leaves (p < 0.05, Table S2.1). As a consequence, the colder the winters and springs, the faster the fire spread (p < 0.05, Table S2.2). However, forest flammability increases significantly with increasing annual, summer or autumn drought (p < 0.05, Table S2.2 in Appendix S2).

Otherwise, partial regressions showed that FLI increases with an increasing proportion of needleleaved trees in the community (r² = 0.17, p < 0.001, Fig. 7a), and that this effect disappeared after accounting for the effect of tree cover (p = 0.12, Fig. 7b). In other words, in the Alps the open canopy structure of needleleaved forests makes them more flammable than the broadleaved.

DISCUSSION

This study proposes the first framework to test how mid-to long-term climate, forest structure and community composition drive the variation in mountain forest understorey flammability, exclusively or in interaction. Interestingly, while precipitation and temperature gradients describe the main axis of variation in community composition (Fig. S1.1), we found that the tree cover explains also a substantial part by leading to the second axis of variation. By controlling light availability in the understorey, the influence of tree cover on low-stature shrub species composition was already observed over the entire Alps [START_REF] Nieto-Lugilde | Tree cover at fine and coarse spatial grains interacts with shade tolerance to shape plant species distributions across the Alps[END_REF]. Further, we found that tree cover and precipitation gradients relate significantly to variations in fire spread and intensity driven by fuel characteristics, thus comforting our predictions. Precisely, this study firstly demonstrates that (H1) understorey flammability varies between mountain ecosystems, mostly according to (H2) the interaction between tree cover, precipitation seasonality and the relative abundance of coniferous trees that jointly determines the variation in fuel structural traits (i.e. fuel load and composition).

Tree cover was found to control the main part of the understorey flammability variance, suggesting that this structural determinant of forests is able to capture most of fuel properties [START_REF] Bradstock | A biogeographic model of fire regimes in Australia: current and future implications[END_REF]. Fire intensity increased with canopy closure until a threshold of tree cover

Author-produced version of the article published in Journal of Biogeography, 2016, 1, 1-12. The original publication is available at http://onlinelibrary.wiley.com/doi/10.1111/jbi.12745/full DOI: 10.1111/jbi.12745 inner continental range where fires are still rare today. Consequently, the flammability drivers we enhanced would benefit from fire prevention and forest management by providing a mechanistic understanding of the areas where the climatic and vegetation features make forests more susceptible to fire damages. In addition this study highlighted the tree cover as a strong leverage from which forest practices may reduce forest flammability in European mountain ecosystems.

Author-produced version of the article published in Journal of Biogeography, 2016, 1, 1-12. The original publication is available at http://onlinelibrary.wiley.com/doi/10.1111/jbi.12745/full DOI: 10.1111/jbi.12745

Author-produced version of the article published in Journal of Biogeography, 2016, 1, 1-12. The original publication is available at http://onlinelibrary.wiley.com/doi/10.1111/jbi.12745/full DOI: 10.1111/jbi.12745

Author-produced version of the article published in Journal of Biogeography, 2016, 1, 1-12. The original publication is available at http://onlinelibrary.wiley.com/doi/10.1111/jbi.12745/full DOI: 10.1111/jbi.12745

Author-produced version of the article published in Journal of Biogeography, 2016, 1, 1-12. The original publication is available at http://onlinelibrary.wiley.com/doi/10.1111/jbi.12745/full DOI: 10.1111/jbi.12745

Supporting Information

Additional Supporting Information may be found in the online version of this article: Appendix S1: Classifying mountain forest ecosystems of the French Alps; designing homogeneous fuel models in heterogeneous forests, and quantifying shrub flammable biomass from individual to community level.

Appendix S2: Relative contribution of litter, herb and shrub to total fine fuel load amongst FET; second-order polynomial regressions relating fuel load to tree cover; and Pearson's correlation coefficients between climate and understorey fuel or flammability parameters Appendix S3: Validation exercises of between-ecosystem differences in and tree cover effect on flammability from lightning fire data Biosketches Thibaut Fréjaville hold a Ph.D. in Environmental Sciences from the Aix-Marseille university and the Mediterranean Ecosystems and Risks Research Group of Irstea Aix-en-Provence (France). His research interests centre on the vulnerability to fire of Mediterranean and mountain forests, including the biogeography of community flammability and fire in response to regional changes in climate and human activities.

Thomas Curt is a researcher at the Mediterranean Ecosystems and Risks Research Group of Irstea Aix-en-Provence (France). His main research interests are in fire ecology and in fire risk assessment in Mediterranean areas.

Christopher Carcaillet is professor at the Ecole Pratique des Hautes Etudes, an institute for advanced studies (Paris). His main research focuses on the disturbance chronoecology in relation to climatic changes and plant diversity dynamics in mountain, boreal and Mediterranean areas, testing the function of disturbances on community pattern and biogeochemical cycles, and exploring the role of global changes (climate, social) on perturbation regimes. Author-produced version of the article published in Journal of Biogeography, 2016, 1, 1-12. The original publication is available at http://onlinelibrary.wiley.com/doi/10.1111/jbi.12745/full DOI: 10.1111/jbi.12745