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 non-linear 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 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.