Modelling interannual variation in the spring and autumn land surface phenology of the European forest

This research reveals new insights into the weather drivers of interannual variation in land surface phenology (LSP) across the entire European forest, while at the same time establishes a new conceptual framework for predictive modelling of LSP. Specifically, the random-forest (RF) method, a multivariate, spatially non-stationary and nonlinear machine learning approach, was introduced for phenological modelling across very large areas and across multiple years simultaneously: the typical case for satellite-observed LSP. The RF model was fitted to the relation between LSP interannual variation and numerous climate predictor variables computed at biologically relevant rather than humanimposed temporal scales. In addition, the legacy effect of an advanced or delayed spring on autumn phenology was explored. The RF models explained 81 and 62 % of the variance in the spring and autumn LSP interannual variation, with relative errors of 10 and 20 %, respectively: a level of precision that has until now been unobtainable at the continental scale. Multivariate linear regression models explained only 36 and 25 %, respectively. It also allowed identification of the main drivers of the interannual variation in LSP through its estimation of variable importance. This research, thus, shows an alternative to the hitherto applied linear regression approaches for modelling LSP and paves the way for further scientific investigation based on machine learning methods.

introduced for phenological modelling across very large areas and across multiple years simultaneously: the typical case for satellite-observed LSP.The RF model was fitted to the relation between LSP interannual variation and numerous climate predictor variables computed at biologically-relevant rather than human-imposed temporal scales.In addition, the legacy effect of an advanced or delayed spring on autumn phenology was explored.The RF models explained 81% and 62% of the variance in the spring and autumn LSP interannual variation, with relative errors of 10% and 20%, respectively: a level of precision that has until now been unobtainable at the continental scale.Multivariate linear regression models explained only 36% and 25%, respectively.It also allowed identification of the main drivers of the interannual variation in LSP through its estimation of variable importance.This research, thus, shows an alternative to the hitherto applied linear regression approaches for modelling LSP and paves the way for further scientific investigation based on machine learning methods.

Introduction
Vegetation phenology has emerged as an important focus for scientific research in the last few decades.The interest in vegetation phenology is twofold: inter-annual recording of the timing of phenological events allows quantification of the impacts of climate change on vegetation; and a greater understanding of phenological responses enables meaningful projections of how ecosystems will respond to future changes in climate (Menzel, 2002;Morisette et al., 2008;Peñuelas, 2009;Peñuelas and Filella, 2001).Although different approaches have been devised for the study of vegetation phenology (Rafferty et al., 2013), the characterisation and modelling of vegetation phenology at global or regional scales has been undertaken mainly through the use of long-term time-series of satellite-sensor vegetation indices (termed land surface phenology, LSP, to reflect that satellite-observed phenology includes all land covers).Most studies of LSP analyse trends in phenological events across years (Delbart et al., 2008;Jeganathan et al., 2014;Jeong et al., 2011;Karlsen et al., 2007;Myneni et al., 1997), but more recent studies present process-based models to uncover cause-effect relationships between long-term trends in phenology and its key driving variables (Ivits et al., 012;Maignan et al., 2008a;Maignan et al., 2008b;Stöckli et al., 2011;Stöckli et al., 2008;Yu et al., 2015;Zhou et al., 2001).This last group of studies focuses on trends in phenology produced by trends in weather (mainly warming).However, interannual variation in LSP arising as a consequence of the inter-annual variability in weather are less studied (Cook et al., 2005;De Beurs and Henebry, 2008;Menzel et al., 2005;Post and Stenseth, 1999;Zhang et al., 2004), with modelbased studies of this phenomenon being scarce (van Vliet, 2010).
A higher frequency in the occurrence of extreme weather events has been observed in Europe, especially for summer temperatures (Barriopedro et al., 2011;Luterbacher et al., 2004).The summers of 2003 and 2010 in western and eastern Europe, respectively, were the warmest in the last 500 years (Barriopedro et al., 2011).Species and ecosystems respond more rapidly to these anomalies in weather than average climatic changes in most climatic scenarios (Zhao et al., 2013).Maignan et al. (2008b) and Rutishauser et al. (2008) reported that the LSP greening occurred 10 days earlier in 2007 than the average over the past three decades as a consequence of an exceptionally mild winter and spring.The study of the impacts of extreme inter-annual weather events on vegetation through the modelling of interannual variation in spring and autumn phenologies can increase our knowledge about climate-driven changes in phenology, acting as natural experiments in climate change scenarios (Rafferty et al., 2013).On the other hand, the modelling of LSP has been less explored compared to the modelling of individual plant species, and there are many aspects that remain to be understood, which limits comprehensive understanding of LSP and, therefore, of phenology at regional or global scales.
A more complete modelling of LSP considering the inter-annual variation across large areas would include the capacity to interpret observations and make meaningful projections in relation to disturbances and their subsequent impacts (Morisette et al., 2008).
Modelling efforts to characterize LSP have generally relied on functions (usually linear) of meteorological drivers, such as average temperature and precipitation (Ivits et al., 2012), growing degree days (GDD) (de Beurs and Henebry, 2005), light and temperature (Stöckli et al., 2011), minimum temperature, photoperiod, vapour pressure deficit (Jolly et al., 2005;Stöckli et al., 2008), or minimum relative humidity (Brown and de Beurs, 2008).However, there is lack of understanding on number of important aspects, such us the multivariate influence of meteorological variables (temperature, precipitation, solar radiation) driving phenology, or the effect of additional drivers in the modelling of autumnal phenophases (Morisette et al., 2008).For instance, Fu et al. (2014) found a "cause-effect relationship" between an earlier leaf senescence and an earlier spring flushing in leaves of warmed samples of Fagus sylvatica and Quercus robur.This legacy effect of spring phenology has been reported in recent studies using modified environments and plant species, but it has not been studied using LSP data.This latter aspect is particularly pertinent for studies that focus on interannual variation in phenology and could potentially contribute to increased knowledge of how climate change is affecting autumn phenology.On the other hand, many studies investigating the sensitivity of phenological events to climate variation use calendar seasonal or monthly mean climatic variables, which operate on fixed human calendar scales with a start date of 1 st of January (Maignan et al., 2008b), instead of using biological scales, for example, time relative to the growing phase of plants (Pau et al., 2011).However, the modelling of interannual variation in LSP considering its potentially complicated relationship with climate in a multidimensional feature space (i.e.high number of multivariate weather drivers) might not be possible using traditional linear regression models (de Beurs and Henebry, 2005).In this sense, phenological modelling may benefit from machine learning techniques such as the Random Forest (RF) method (Breiman, 2001), reducing uncertainties and bias (Zhao et al., 2013).RFs have the potential to identify and model the complex non-linear relationships between phenology and climate, being able to handle a large number of predictors and determine their importance in explaining phenology.RFs has been applied with very promising results to other fields of ecology and biological sciences (Archibald et al., 2009;Darling et al., 2012;Lawler et al., 2006), as well as to the simulation of phenological shifts under different climatic change scenarios (Lebourgeois et al., 2010), but the potential for modelling climate-driven interannual variation in phenology is still to be explored.
Understanding the effect of inter-annual weather variation on LSP is an essential step to establish a plausible link between recent climate variability and vegetation phenological responses at global or regional scales, and importantly to make reliable forecasts about future vegetation responses to different future climatic scenarios.The aim of this study is, therefore, to provide an explanation of the observed interannual variation in LSP of the entire European forest during the last decade, identifying the main weather drivers for spring and autumn at the continental scale.Our research offers new insights into the study of LSP by modelling the climate-driven past interannual variation in phenology, rather than trends, and using innovative multivariate non-linear machine learning techniques to evaluate multiple weather predictors at biological scales, and non-weather predictors such as the legacy effect of the date of spring onset in leaf senescence.Climate predictors used range from 30 days average values of temperature variables (max, min and avg) such as precipitation, short wave radiation and day length; trimestral cumulated values such as growing degree days or chilling requirements, among others; to the date of specific events such as the first freeze or the last freeze.Moreover, we considered flexible biological time scales in the analysis between weather and phenological events rather than calendar months.

Data
Three sources of data were used for this research: i) Satellite sensor derived temporal composites of MERIS Terrestrial Chlorophyll Index (MTCI), ii) temperature and precipitation data from the European Climate Assessment and Data (ECA&D) project (http: //www.ecad.eu)and iii) surface radiation daylight (DAL; w/m 2 ) data and surface incoming shortwave (SIS; w/m 2 ) radiation data from the Climate Monitoring Satellite Application Facilities (CM SAF, http://www.cmsaf.eu).
We used weekly composites of MTCI data at 1 km spatial resolution from 2002 to 2012.This dataset was supplied by the European Space Agency and processed by Airbus Defence and Space.Daily temperature (mean, minimum and maximum) and daily precipitation data were derived from the European Climate Assessment & Dataset (ECA&D) time-series (version 10.0) with spatial resolution of 0.25° ×0.25°, covering the period from 2002 to 2011 (Haylock et al., 2008).The CM SAF DAL version CDR v001 (Müller and Trentmann, 2013) and SIS version CDR v002 (Posselt et al., 2012;Posselt et al., 2011) were derived from Meteosat satellite sensors at a spatial resolution of 0.05° x0.05° covering the same period as ECA&D.

Phenology extraction and interannual variation in LSP computation
The time-series of MERIS MTCI data was used to estimate both the onset of greenness (OG) and end of senescence (EOS) from 2003 to 2011.Data for every estimation year considered 1.5 years of data (from October in the previous year to July in the next year) because the annual pattern of vegetation growth in some parts of Europe spans across calendar years and, hence, insufficient information about LSP is captured using a single year of data.The yearly values of OG and EOS were estimated for each image pixel of the study area using the methodology described in Dash et al. (2010).This methodology consists of two major procedures: data smoothing and LSP estimation (Figure 2a).
Smoothed MTCI time-series data were obtained using a discrete Fourier transform because of its advantage of requiring fewer user-defined parameters compared to other methods (Atkinson et al., 2012).The peak in the annual profile was defined as a point on the phenological curve where the first derivative changes sign from positive to negative.Next, the derived data were searched backward and forward departing from the maximum annual peak to estimate the OG and EOS, respectively.OG was defined as a valley at the beginning of the growing season point (a change in derivative value from positive to negative) and EOS was defined as a valley point occurring at the decaying end of a phenology cycle (a change in derivative value from negative to positive).These satellitederived LSP estimates were compared to ground observations of the thousands of deciduous tree phenology records of the Pan European Phenology network (PEP725) (Rodriguez-Galiano et al., 2015a).This comparison resulted in a large spatio-temporal correlation of the phenology estimates with the spring phenophase (OG vs leaf unfolding; pseudo-R 2 =0.70) and autumn phenophase (EOS vs autumnal colouring; pseudo-R 2 =0.71).
Z-score values during the study period were used as a proxy to measure interannual variation in the LSP parameters.The z-score values for a given year were defined as the difference from the multi-year mean, normalized by the standard deviation across years.
The value of the targeted year was excluded in the computation of multiyear mean to enhance the inter-annual variation (Saleska et al., 2007).The spatio-temporal distribution of spring and autumn LSP z-score values is shown in Figures S1 and S2 of the supporting information, respectively.
To match the spatial resolution of the ECA&D dataset, the LSP z-score values for each year were resampled to a spatial resolution of 0.25°×0.25°by calculating the median of all the LSP z-score values within this area after excluding the areas with fewer than 50 LSP estimates and the non-forest pixels according to the Globcover2005 and Globcover2009 land cover maps (http://due.esrin.esa.int/globcover/).Only LSP estimates with complete temporal coverage (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) were included in the analysis to reduce the likelihood of natural and human disturbances (Potter et al., 2003).Globcover was selected for its greater consistency with the MERIS MTCI time-series and its high geolocational accuracy (<150 m) (Bicheron et al., 2011).

Computation of weather predictors
A suite of weather predictors were computed for each 0.25 ×0.25° grid cell associated with the occurrence of positive or negative z-score values in LSP based on the ECA&D and CM SAF datasets (see Table 1).The predictors include temporal average values of temperature variables (Tmax, Tmin and Tavg), precipitation, DAL and SIS; temporal cumulated predictors such as growing degree days, chilling, precipitation, SIS and DAL; and the date of specific events such as the onset of greenness (legacy effect for autumn phenology modelling) the first freeze or the last freeze, as well as the difference between both dates (freeze period) for the modelling of autumn only.Growing degree days were computed using temperature thresholds of 0° and 5°.
Chilling requirements were computed as the sum of negative temperatures (temperatures below 0°).Freeze was defined as dates with minimum temperatures lower than -2° (Schwartz et al., 2006).
The different weather predictors were computed based on the 30 and 90 days previous to the day of the year (DOY) of the z-score values in OG and EOS (Figure 2b) following Schwartz et al. (2006) and Menzel et al. (2006), who found that most phenophases of plant observations in Europe correlated significantly with weather predictors representing the month of onset and the two preceding months.The chilling requirements for spring modelling and freeze predictors were an exception, as the period for its computation starts 90 days prior to the OG.Relative differences between each predictor and its multi-year average for the same period were computed to capture the inter-annual variability in climate variables at the pixel level for every predictor and to facilitate the modelling of climate-driven variation in phenology (Table 1).

Modelling interannual variation in LSP
Conventional statistical models such as linear regression might be inappropriate for investigating the drivers of interannual variation in phenology because many of the relationships are likely to be non-linear (De Beurs and Henebry, 2008).In this sense, machine learning methods have emerged as complementary alternatives to conventional statistical techniques.Within the branch of machine learning techniques, regression trees are particularly suitable when compared to global single predictive models, allowing for multiple regression models using recursive partitioning (Breiman, 1984).Assembling a single global model might not be representative of LSP of the entire European continent, when there are many climatic drivers which interact in complicated, non-linear ways and may vary spatially and temporally.
For the purpose of this paper, an alternative approach is to sub-divide, or partition, the data space into more homogeneous regions of similar climates and ecological factors.
Regression trees use a sum of squares criterion to split the data into successively more homogeneous subsets contained at many different structural units called nodes.Each of the terminal nodes, has attached to it a simple regression which applies in that node only.Therefore, different regressions can be fitted to different data subsets within one single regression tree, which can represent different responses controlled by different drivers (Archibald et al., 2009;Lawler et al., 2006).Additionally, the performance of multiple regression trees can be combined to increase the predictive ability of a single regression tree model, following the Random Forest technique (Figure 3).The RF method is an innovative machine learning approach that can perform multivariate non-linear regression, combining the performance of numerous regression tree algorithms to predict the interannual variation in OG and EOS.More The Random Forest method was applied to phenological modelling across very large areas and across multiple years simultaneously: the typical case for satellite-observed LSP.The RF model was fitted to the relation between LSP interannual variation and numerous climate predictor variables computed at biologically-relevant rather than human-imposed temporal scales.We restricted our climate data choices to daily data (average, minimum and maximum temperatures, precipitation and radiation) to account for integrative forcing (that is, growing degree days, chilling requirements as well as cumulative precipitation and radiation), computed from the exact day of the phenological event backwards, rather than using the calendar months.
The locations with z-score in LSP greater than 1 (positive and negative) were selected to build a RF predictive model on OG and EOS.Z-score values of OG or EOS for each year were combined together with the different weather predictors.The z-score values in OG were assessed as an extra predictor to evaluate the legacy effect of an advanced or delayed spring in the modelling of EOS.The values of these variables at the selected years and locations (spatiotemporal model) were combined into a set of input feature vectors (3900 feature vectors for the spring model and 3124 for autumn) as an input to the RF algorithm.These feature vectors were divided equally into two subsets, one for the training of the models (inbag) and one as an additional test to the one internally computed by RF (out of bag; oob) to evaluate performance.RF models composed of 2000 trees were grown using different subsets of predictors, varying the number of random predictors from 1 to 9. The Random Forest method within the package implemented in the R statistical software was used to build the different models (Liaw and Wiener, 2002).

Selection of the most important predictors
The RF method can use the oob subset to estimate the relative importance of each predictor in the model.This property is especially useful for the present research, but also for other multivariate biological studies, where it is important to know the physical drivers of the phenomenon under investigation (Archibald et al., 2009;Lawler et al., 2006).However, the inclusion of different measures of weather predictors may imply a large increase in the dimensionality of the datasets being used, as these variables are obtained by applying multiple functions or measures to the temperature, precipitation and radiation time-series.On the one hand, more information may be useful for the modelling process; on the other hand, an excessive number of correlated predictors or features can overwhelm the expected increase in accuracy and may introduce additional complexity limiting the ability of the method to point to possible cause-effect relationships between interannual variation in phenology and their drivers, making interpretation challenging.
A feature selection approach, based on the ability of the RF to assess the relative importance of the predictors, was used to identify the minimum number of drivers which can better explain spring or autumn interannual variation in phenology.To assess the importance of each weather predictor, the RF switches one of the input predictors while keeping the rest constant, and it reevaluates the performance of the model measuring the decrease in node impurity (Breiman, 2001).The differences were averaged over all 2000 trees to compute the general drivers for the interannual variation in Europe.However, different subsets of variables could be used to characterize different climates and ecological factors at every single regression tree model or node (see previous section).In order to reduce the number of drivers the least important predictor was removed iteratively at different steps.Then, a 5-fold cross-validation was applied to obtain a stable estimate of the error of the model built after predictor deletions.Finally, the model with a better trade-off between number of predictors and error was chosen as the basis for interpreting the likely drivers of interannual variation in phenology.

Results
Numerous models were built on the basis of different predictor combinations considering different temporal windows prior to the spring and autumn phenological events (see section "computation of weather predictors").The percentage of variation (pseudo-R 2 ) explained by different weather-LSP models is shown in the supplementary information (Table S1, S2 and    S3).No previous studies have investigated in depth the parametrization of GDD for LSP and climate inter-comparison, unlike for ground phenological studies (Snyder et al., 1999).
Although, we did not carry out an exhaustive analysis of the optimum GDD parametrization, our results showed a systematic pattern in spring models, presenting slightly larger pseudo-R 2 for models which used 0⁰ C as a threshold for the computation of GDD (rather than 5⁰ C).
Regarding, the length of the temporal windows for weather function computation, spring models using 30 and 90 days for the computation of averaged and cumulative functions were more accurate, whereas for autumn models with 90 day-averaged predictors outperformed the rest.
The main drivers of interannual variation in LSP were identified through the application of a feature selection procedure (see section "selection of the most important predictors").Spring models were more accurate than autumn, with median relative error values of 10% to 27% (12 to 1 predictor), versus 26% to 60% of autumn (14 to 1 predictor).Figure 4 shows the pseudo-R 2 of the models as well as the relative importance of each predictor.Spring models (explained a percentage of the variance up to 81% (Figure 4a), whereas autumn explained up to 61% (Figure 4b).Cook et al. (2005), using a modelled based on GDD only, explained 63% on the variance of onset date for mixed and boreal forest.Figure 5 shows the relative error in the prediction of different models after removing the least important predictor.Regarding the relative importance of the drivers, the same ranking in importance was observed within the different models of each phenophase, which reflected the stability in the RF importance estimation, and a high reliability of the results (Figure 4).To interpret the main weather drivers of the interannual variation in phenology, simplified models with reduced number of predictors were selected for spring and autumn (see section 3.5), respectively.The spring model was composed of 6 predictors (pseudo-R 2 =0.77 and median relative error of 10%) and the autumn model of 5 predictors (pseudo-R 2 =0.59 and median relative error of 28%) (Figure 6).Our results suggest that interannual variation in the onset on greenness (LSP) of temperate forest species are driven mainly by the daily temperature of the 30 days prior to onset (but not necessarily the GDD), with the most important driver being the minimum temperature.
Photoperiod was also important, the most accurate empirical prediction was obtained by a combined temperature-radiation forcing, integrating the SIS of the previous 90 days.For senescence, temperature was suggested to be more important than photoperiod in controlling the senescence process (Archetti et al., 2013;Jeong and Medvigy, 2014;Vitasse et al., 2009;Yang et al., 2012), with the most important drivers being the date of the first freeze and the accumulation of chilling temperatures.However, we did not observe a legacy effect of a much earlier or later spring onset on the date of senescence.Autumn models that included the interannual variation (z-score values) in the onset of greenness did not outperform the remaining models (see Table S2 and S3 in supplementary information) and the relative importance was low in comparison with other drivers.

Discussion
The selection and computation of the weather predictors is an important step of phenological modelling.Most of studies on the sensitivity of phenological events to climate used human calendar scales, that is, seasonal or monthly calendar mean or cumulative climate predictors (Maignan et al., 2008a;Maignan et al., 2008b;Menzel et al., 2006;Schwartz et al., 2006), overlooking the importance of biological time-scales in phenology.However, with the increased availability of daily weather datasets, current and future studies might benefit from the use of daily information to model the drivers of plants' circadian time-scales (Pau et al., 2011).Our study advanced the modelling of vegetation phenology by improving the temporal matching between LSP interannual variation and the preceding weather conditions by analysing daily data at biological scales.Regarding, the length of the temporal windows for weather function computation, Menzel et al. (2006) showed that most phenological phases of plant species in Europe correlate significantly with mean temperatures of the month of onset and the two preceding months.However, in our study, when end of senescence was considered, a consistent divergent effect was observed between spring and autumn.Autumn phenophases might be driven by longer-term changes in weather, while for spring the average conditions of the 30 days previous to the date of onset play a more important role (Table S1, S2 and S3 in supplementary information).From a computational point of view, considering larger temporal windows for calculating averages would induce a smoothing effect, degrading the information in the predictors, whereas cumulative functions such as GDD or chilling requirements would not be affected by this effect.However, we observed a divergent response between spring and autumn and consistent throughout the models of each phenophase suggests that a biological explanation for this phenomenon might be plausible.
Understanding the drivers of interannual variation in LSP amidst background inter-annual variation is a critical aspect of global change science (de Beurs and Henebry, 2005;Zhao et al., 2013).To this end, the RF method is particularly pertinent, as it allows the assessment of the importance of the predictors (Figure 4).Our findings reveal that the accuracy of growing degree day-based models might be overestimated using linear regression models and that non-linear multivariate relationships between temperature (especially minimum temperature) and radiation are needed to describe the relations between phenology and weather drivers.This supports the findings of Stöckli et al. (2011) who explained temperate phenology using a combination of light and temperature.The highlighted importance of minimum temperatures might be related to the fact that minimum temperature is a better indicator of weather changes than either the average or maximum temperature (Duncan et al., 2014;Jolly et al., 2005).
Regarding GDD, although it has been applied extensively to predict vegetation phenophases , it is currently debated whether such models can detect when multiple environmental drivers are required to initiate a phenological event, or detect drivers that are relatively static across time, such as photoperiod (Stöckli et al. 2011).Our results reveal that multiple environmental drivers are required to initiate phenological events of Europe and also showed that the role of GDD alone in driving spring phenology might be overestimated due to an over-reliance on linear models.GDD had the largest linear association with vegetation phenology interannual variation, while the linear correlation between LSP and others drivers that were revealed as very important by the RF was small (see Tables 1 and 2).A simple linear analysis between GDD and phenology could ignore complex non-linear associations between phenology and predictors as well as synergies between weather drivers.Regarding the senescence phase, the autumn models had a weaker predictive power compared the spring models.There is still lack of clear understanding of mechanism autumn senescence, however, temperature, and particularly the dates of freeze, has been suggested as major driver for autumn phenology.
The RF method provided an important alternative over simple, but less accurate analysis based on linear regression for the analysis of interannual variation in spring and autumn phenology.
A further comparison with a linear regression analysis suggested that there might be a nonlinear relationship between the interannual variation in LSP and the weather drivers.
Multivariate linear regression models were also fitted from the same combination of predictors selected as optimal by Random Forest.Multivariate linear models explained only 36% and 26% of the variance in spring and autumn phenology interannual variation across the continental scale.Additionally, a linear regression between predicted values from RF and observed interannual variation in phenology produced R 2 values equal to 0.90 and 0.68 for spring and autumn LSP interannual variation, respectively (Figure 6a and 6b).On the other hand, the correlations between the predictions of linear regression models and observations were much weaker, with R 2 values of 0.39 and 0.25 (Figure 6c and 6d).Linear models under-predicted a delay in the phenophases (positive z-score values) and over-predicted the advances (negative z-score values).The spatial distribution of the relative errors for RF and multivariate linear regression is shown in Figures S3 to S6 of the supporting information.The relative errors of the latter were significantly higher.Additionally, the residuals seemed not to be homoscedastic suggesting that linear models might not be able to deal with the complex patterns between LSP and climate patterns at multiple locations and times, integrating them into a unique overall model.
A new approach to model interannual variation in LSP was presented in this paper based on the application of the RF model to a set of climate predictors at biological scales.This new modelling technique has numerous advantages for the modelling of climate-driven interannual variation in LSP.It is a non-parametric multivariate method which allows for non-linear relationships between (compared to traditional linear models) phenology and climate and can consider a large number of weather predictors in the modelling process.This provides potential opportunity to capture the impact of all possible environmental/weather drivers on vegetation phenology.The proposed method can recognize complex patterns between LSP and climate at multiple locations and times, integrating them into a unique overall model, rather than generating multiple models over a geographical area and for different years.Additionally it is data-driven, which means that there is no need to incorporate previous knowledge about the specific responses of vegetation to different predominant weather controls (i.e.temperature, rainfall, and photoperiod), allowing weather drivers to automatically shift both temporally and spatially.Therefore, it is highly generalizable, being applicable to different biogeographical regions where the phenology is controlled by different factors.This flexibility or generalization capacity of RF models to transition from one driver to another without the need for a model change also promotes its application to different climate change scenarios.We succeeded in modelling the interannual variation in LSP phenology as observed from satellite-sensors in the European Forest, while using the same type of input data, the same model, and the same model parameters for the entire European continent.
Table 1.Predictors used in the modelling of the interannual variation in LSP.* predicted over a period of 90 days.** predicted over a period of the 30 and 90 days previous to the date of the z-score value.
details regarding the performance and the specific characteristics of a RF model can be seen in Rodriguez-Galiano et al. (2015b); Rodriguez-Galiano et al. (2014), and Figure 3.

Figure 2 .
Figure 2. Flow-chart illustrating the methodology.A) Phenology extraction and interannual variation in LSP computation.B) Computation of weather predictors.C) Modelling of interannual variation in phenology.

Figure 3 .
Figure 3.The flowchart of Random Forest for regression (adapted from Rodriguez-Galiano et al. 2015b).The RF method receives a subset of input vectors (n), made up of one phenology zscore value and the values of the corresponding weather predictors for a given location and year.RF builds a number K of regression trees making them grow from different training datasubsets, resampling randomly the original dataset with replacement.Hence, most data will be used multiple times in different models.On the other hand, when the RF makes a tree grow, it uses the best predictor within a subset of predictors (m) which has been selected randomly from the overall set of input predictors.These especial characteristics of RF confer a greater prediction stability and accuracy and, at the same time, avoid the correlation of the different RTs, increasing the diversity of patterns that can be learnt from data.The multiple predictions of all k RTs for a given vector used as training are then averaged to obtain a unique estimation of the phenology z-score value.

Figure 4 .
Figure 4. Relative importance of each independent variable in predicting phenology interannual variation in Europe.Different models derived from the feature selection approach are represented in each column.Numbers given over each column represent the coefficient determination of each model.Plots at the top and bottom represent the spring (a) and autumn interannual variation in LSP (b), respectively.The names of predictors follows the notation: Prefix M and C represent the mean and cumulated functions; TX, TN and TG: maximum, minimum and average temperature, respectively; PP: precipitation; SIS: surface incoming shortwave radiation; DAL: surface radiation daylight; GDD: growing degree days; CHIL: chilling requirements; FF, LF and PF: first, last and period of freeze, respectively.

Figure 5 .
Figure 5. Relative error of the models fitted as a result of the feature selection approach.Median (interior horizontal line), mean (interior square), 1% and 99% quantiles (edge of boxes), range (extremes).Relative errors were calculated for the prediction of 1,974 and 1,576 independent observations for spring (a) and autumn (b), respectively.See previous figure for the weather predictor variables in the models, as shown in the x-axis.

Figure 6 .
Figure 6.Scatterplots between observed anomalies in LSP and the predictions calculated using a selection of weather predictors (see Figure 2 and Figure 3).Plots for spring phenology are shown on the left panel (blue; a, c) and autumn on the right (red; b, d).Random Forest predictions are given in the upper panel (a, b) and those of the linear regression in the bottom (c, d) panel.The dashed lines represent an exact 1:1 relationship (expected fitting), the solid lines show a linear regression of these data.The explained variances (percentage R 2 ) and RMSE values are 90% and 0.43 (spring Random Forest model), 68% and 0.92 (autumn Random Forest model), 39% and 1.04 (spring Linear model) and 25% 1.40 (autumn linear model).

Table 2 .
Correlations between the predictors used in the modelling of spring interannual variation in LSP.Significant correlations between the 1 anomalies and the predictors are given in bold (p < 0.05).

Table 3 .
Correlations between the predictors used in the modelling of autumn interannual variation in LSP.Significant correlations between the 1 anomalies and the predictors are given in bold (p < 0.05).