Patterns and trends in the spatial heterogeneity of land surface phenology of global forests

Land-surface phenology is a widely used indicator of how terrestrial ecosystems respond to environmental change. The spatial variability of this plant functional trait has also been advocated as an indicator of the functional composition of ecosystems. However, a global-scale assessment of spatial patterns in the spatial heterogeneity of forest phenology is currently lacking. To address this knowledge gap, we developed an index based on satellite retrievals and used it to quantify phenological diversity across global forest biomes. We show that there is considerable variation in phenological diversity among biomes, with the highest overall levels occurring in arid and temperate regions. An analysis of the drivers of the spatial patterns revealed that temperature-related factors primarily determine the variation in phenological diversity. Notably, temperature seasonality and mean annual temperature emerged as the most significant variables in explaining this global-scale variability. Furthermore, an assessment of temporal changes over an 18-year period revealed strong climate-driven shifts of phenological diversity in boreal and arid regions, suggesting that there may be an ongoing widespread homogenisation of land surface phenology within forest ecosystems. Our findings ultimately contribute to the development of a novel Essential Biodiversity Variable, which may enable scientists and practitioners to quantify the functional composition of ecosystems at unprecedented spatial and temporal scales.


Introduction
Forest ecosystems provide a range of functions, which are vital for supporting the provision of goods and services upon which human societies depend.To understand how they may respond to future environmental change, it is of utmost importance to quantify the spatial and temporal distribution of their functional structure.To date, functional biogeography (Reichstein et al 2014, Violle et al 2014) has focused on a wide range of issues, including the analysis of temporal changes in a number of functional properties (Musavi et al 2017, Soranno et al 2019, Vaglio Laurin et al 2020).However, the study of spatial heterogeneity in functional properties as a means of understanding their functional composition has received much less attention (Jetz et al 2016).Quantifying spatial heterogeneity in functional properties has the potential to further improve our understanding of the role of ecosystem functional composition in modulating responses to environmental change.Additionally, it may shed light on the mechanisms driving ecosystem performance and resilience across large geographical gradients.(Walter et al 2017, Dronova andTaddeo 2022).
Phenology is a key plant trait and a functional property of terrestrial ecosystems, which can be used to describe the seasonal dynamics of vegetated land.While phenology is generally measurable at individual and species-level, it can also be measured as an aggregate property for ecosystems.Land surface phenology is typically estimated from data acquired using space-borne optical sensors (Helman 2018).This offers a complementary metric, compared to that provided using ground-based data, and is capable of capturing multiple phenological trajectories for a given area or different vegetation types.Land surface phenology is a key indicator of global change that both responds to and influences weather and climate (Richardson et al 2013) and has been extensively used to study long-term vegetation dynamics at a variety of spatial and temporal scales (Buitenwerf et al 2015, Forkel et al 2015, Garonna et al 2016).
This has been largely achieved through the analysis of time series of vegetation indices, which can be used to quantify inter-annual changes in the timing and intensity of vegetation activity, and from which a variety of phenological metrics can be derived (Helman 2018).However, more recently, the use of phenological traits has been advocated as a promising tool for quantifying plant diversity at large spatial scales (Viña et al 2016, Dronova andTaddeo 2022).Furthermore, phenological traits have been proposed as an important class of Essential Biodiversity Variables (EBV).However, an assessment of the spatial and temporal variability phenological diversity across broad ecological gradients is currently lacking.Phenological diversity, quantified using land surface phenology data, may reflect functional aspects of ecological communities related to resource utilisation and niche differentiation in time (Dronova and Taddeo 2022).Indices describing spatial heterogeneity in land surface phenology have therefore the potential to be informative about the functional composition of ecosystems and their changes.Satellite-based phenodiversity metrics may thus be complementary to existing heterogeneity metrics based on snapshot observations from spectral indices or raw reflectance data from spaceborne passive sensors (Maeda et al 2014, Dronova and Taddeo 2022, Fassnacht et al 2022).
In this study, we quantify global patterns in forest phenological diversity across space and time using longterm satellite remote sensing data.To do so, we develop a novel index that quantifies the spatial variability in Land Surface Phenology for global forested areas, based on cumulated Normalized Difference Vegetation Index (NDVI) time-series data.Specifically, phenological diversity is computed as the spatial variance of the phenology of individual phenological trajectories for a given area (see Methods).Our index is designed to minimise the impact of noise in the observation records due to the uncertainty of satellite retrievals (i.e.generated by aerosol load, cloud cover, and co-registration errors).The information produced is temporally explicit, enabling us to draw a comprehensive picture of interannual changes and trends in phenological diversity.The specific objectives of this study are twofold: (i) to quantify the primary drivers of spatial patterns in global forest phenological diversity and (ii) to quantify changes in phenological diversity that have occurred during the period 2003-2020 while also quantifying their sensitivity to two key climatic drivers.

Remote sensing data
The remote sensing data used in this study consists of MODIS multi-day composites surface reflectance data with a spatial resolution of 500m and a temporal resolution of 8 days (Justice et al 2002).We decided to use MODIS data because of the consistency of data derived from a single sensor, ensuring uniformity and reliability in our time-series analysis.This consistency is essential for accurately detecting and interpreting ecological changes over time, minimizing the variations and uncertainties often associated with data from multiple sensors.While the 18-year span of our dataset may seem relatively short, it aligns well with similar studies in the field and is sufficiently long to capture significant ecological changes, particularly in the context of recent accelerated climate trends (Moon et al 2021, Gao et al 2023).The reflectance dataset needs to be combined with the quality state flags available at ~1 km to remove undesired pixels.The corresponding flags are used to remove pixels that are suspected to contain clouds and cloud shadows, or if the pixels are adjacent to suspected clouds.Pixels located in oceans beyond ocean coastlines or shallow oceans are also removed.This was achieved by using the StateQA band from the MODIS dataset to filter out pixels affected by clouds, cloud shadows, adjacency to clouds, as well as those located in oceans beyond coastlines or in shallow ocean areas.Additionally, the filtering process used internal flags from the StateQA to exclude pixels containing fire and snow/ice.The input NDVI data was calculated using bands 2 (near-infrared) and 1 (red).Because gaps within the time series may have a strong impact the estimation of phenology, we used a gap-filling procedure.In brief, gaps at the 8-day level were filled by means of linear interpolation, using a window of 8 observations.Remaining gaps due to persistent clouds were gap-filled using the 5 -year mean values of the NDVI (figure S1, figure S6).

Phenological diversity metric
Prior to the calculations of our phenological diversity metric, we masked out non-forested areas and areas of forest that had changed into other land cover types over the period of 2000-2018.This mask was created using the ESA land cover cci product (http://maps.elie.ucl.ac.be/CCI/viewer/download/ESACCI-LC-Ph2-PUGv2_2.0.pdf see also figure S1) using annual land cover data, ranging from 2000-2018.We chose the ESACCI land cover product to identify areas where forests had been stable through time over other land cover products because of its usefulness in climate change studies.The ESA CCI land cover dataset, produced as part of the ESA Climate Change Initiative, focuses on generating consistent global land cover maps.The land cover product offers a stable and reliable dataset, which is crucial for identifying areas where forests remain stable over time.
Phenology is typically characterised by a series of metrics that describe specific features of the phenological cycle, such as the start of the season, end of the season, and growing season length.However, the use of these classical metrics to define a global phenological diversity index is limited by several factors including (1) their limited applicability to areas with multiple growing seasons (e.g., meadows or semi-arid regions) or continuous growth (e.g., humid tropical forests or boreal forests) (2) the subjectivity involved in combining various metrics into a single index (3) significant uncertainty in computing the metrics in regions with persistent cloud cover during phenological transitions, requiring the need for temporal interpolation of time series (e.g., using spline or double logistic functions) and (4) the high noise-to signal ratio in heterogeneous areas, where satellite retrievals are affected by a higher random error due to spatial sampling of the surface (e.g., particularly relevant for multiangular sensors with varying footprints like MODIS).
To overcome the limitations of classical phenological metrics, for an assessment of phenological diversity, we developed and applied a novel metric based on the following key features: (1) The phenology of each pixel is characterised by the shape of the cumulative time series of NDVI.The trajectory of the cumulative NDVI has the advantage to (a) reduce the importance of the random error on the metrics extraction due to the retrieval algorithm, the atmospheric conditions and the spatial heterogeneity, (b) to condensate in one single dimension all the multiple potential phenological metrics that characterise the annual phenological cycle, and (c) it is a applicable to any biome and phenology strategy (e.g.multiple growing seasons).(2) Normalisation of cumulative NDVI by dividing each element of the time series by the cumulative value at the end of the time series.This normalisation procedure ensures that the cumulative NDVI curves and the derived index only reflect the signal of vegetation phenology of individual pixels and not the spatial variation in NDVI and productivity.
(3) Assessment of the spatial variability of the normalised cumulative NDVI.Once the cumulated and normalised NDVI curves are computed for each MODIS pixel at the resolution of 500 m, the spatial diversity of phenology is assessed within a spatial window with a resolution of ~5 km, equal to 100 pixels.For each time step within a given year (46 time steps, figure 1), the diversity is computed as the deviations of the mean values of the phenological trajectories for all the 500 m × 500 m pixels contained within the 5 km × 5 km spatial window.A graphical description of the algorithm is shown in figure 1.
Below we provide a detailed technical description of the steps involved in the computation of the index: 1. Cumulative sum of the NDVI where CS j i [ ] represents the j-th element of the i-th cumulative sum vector CS i and V j k [ ] represent the j-th element of the k-th vector containing NDVI values for a given year.

NormCS
where CS j k [ ]represents the normalised value of the j-th element of the k-th cumulative sum vector CS k and k w represents the last value of the cumulative sum vector CS k 3. Calculate deviations from the mean value of the cumulative normalised for each time step where t s represents the standard deviation at the t-th time step, d kt represents the deviation of the k-th normalised vector at time step t, t m represents the mean deviation at time step t calculated across the normalised vectors and N represents the total number of normalised vectors.The values obtained from step 3 are then averaged across all the time steps.
Our approach provides a quantitative description of the spatial variability of phenological trajectories for a given area that (1) summarises the diversity of the entire annual phenological cycle in a single metric, (2) works in cloud-persistent areas and (3) is applicable to any ecosystem and phenological strategy.The algorithm was used to produce yearly maps of phenological diversity.All the calculations were performed using Google Earth Engine, a multi-petabyte analysis-ready data catalogue co-located with a high-performance, intrinsically parallel computation service (Gorelick et al 2017).For subsequent analyses, we further filtered the results and retained only ~5km pixels that had a forest cover 0.85.This ensured that only purest forest pixels would be used, while making sure that subsequent analyses would be computationally tractable.

Analysis of spatial drivers of phenological diversity
We first conducted a variance components analysis to explore where most variation in phenological diversity is found.This was achieved by using a linear mixed effect model (LME) with random effect terms for biomes and plant functional type, using the REML (restricted maximum likelihood) criterion .The variance components represent the variance of random terms and the random error of the model.The contribution of each source of variation can quantified via the mixed effect model.
To quantify the drivers of spatial variation in phenological diversity we created a multi-annual average of the yearly global dataset.We analysed the relations between phenological diversity and ten predictors, including variables related to climate, topography and the human pressure on the environment (see table S1).We used long-term average data for five climatic variables from the CHELSA dataset (https://chelsa-climate.org/), including cumulative precipitation, average annual temperature, temperature seasonality (i.e. standard deviation monthly temperatures), precipitation seasonality (i.e.Standard deviation monthly cumulative precipitation.), and aridity.The spatial average and standard deviation (indicating spatial heterogeneity) were calculated from pixel values within each ~5km window previously used for the calculation of phenological diversity.A detailed list of the climatic variables used is available in table S1.We also used a measure of topographic complexity, derived from https://earthenv.org/topography.This measure is calculated using the Shannon Index with a dataset describing the most common geomorphological forms.Additionally, we obtained a measure of human impact from the Human Footprint dataset (https://earthenv.org).This index is a composite metric quantifying human pressure on the environment and is derived from spatially-explicit data for eight human pressures.We used a machine-learning method, Random forests (Breiman 2001), to quantify the relations between the aforementioned predictors (climate, human impact and topography) and phenological diversity (used as dependent variable).We choose the Random Forests algorithm because of its ability to handle non-linear relationships with multidimensional data, as well as its capacity to handle multicollinearity among predictors.Random Forest's ability to handle multicollinearity lies in the random selection of predictor variables at each split in the decision trees.This process is very effective as it crucially allows variables that might otherwise be overshadowed by more dominant, correlated predictors in traditional models to contribute to the ensemble (Cutler et al 2007, Genuer et al 2010).The optimization of the Random Forests model involves tuning a number of hyperparameters.These include the size of the forest (i.e. the number of regression or classification trees) and the randomness or mtry parameter (i.e. the number of variables considered as candidate for growing an individual tree).We specified a fixed number of trees at 500.This decision was made to maximise computational efficiency, as increasing the number of trees beyond a certain threshold typically results in marginal gains in performance at the cost of higher computational demands.The choice of 500 trees balances model accuracy and computational efficiency, ensuring robustness in capturing complex data patterns.To identify the best mtry value, we employed a grid search approach.This method involves systematically testing a range of mtry values and evaluating the model's performance for each.The mtry parameter is crucial in determining the subset of features considered at each split in the trees and therefore significantly influences the diversity and independence of the trees in the forest.By performing a grid search, we aimed to find an mtry value that optimizes model accuracy while preventing overfitting.The optimal mtry value is typically dataset-specific, so this process ensures that our model is finely tuned to the peculiarities of our data.
A 10-fold block cross-validation (Valavi et al 2019) approach was employed to quantify the model's ability to predict phenological diversity.In brief, the dataset was divided in spatially independent blocks, in such a way that all the blocks constituting a fold would be spatially independent.Independent blocks were created using a pre-specified distance.The distance over which observations are independent was determined by constructing the empirical variogram of the residuals of a model initially fitted using the 10 predictors.To evaluate the influence of each predictor in our Random Forests model, we implemented a permutation importance procedure (REF).This method, known for its robustness in feature selection, systematically alters the order of each predictor variable and measures the resulting variation in the model's Mean Squared Error (MSE).Specifically, the permutation process involves randomly shuffling the values of each predictor variable across the dataset while keeping the values of all other variables constant.A significant decrease in R 2 upon permutation indicates a high importance of the variable.This method effectively isolates the impact of each predictor on the model's performance, offering a clear quantification of its relative importance.The process is repeated multiple times for each predictor to ensure statistical robustness.The average decrease in model accuracy across these repetitions was calculated to derive the final importance score for each variable.The response functions of the fitted random forest model was examined using partial dependence plots (PDPs) (Friedman 2001).Twodimensional partial dependence plots were used to identify which predictor-interactions may have an effect on phenological diversity.
To examine local responses of phenological diversity to the predictors we used the LIME (Local agnostic model explanations) methodology (Ribeiro et al 2016).The aim of LIME is to explain how the fitted complex model creates a prediction for a given instance (i.e. a grid cell or other local neighborhood).To this end, for each instance, LIME fits a 'local surrogate' model (a simple model) that approximates the behavior of the complex model for a limited area of the n-dimensional space defined by the predictor variables.We use a ridge regression model for creating LIME surfaces of local coefficients.

Temporal changes in phenological diversity
To analyse changes in phenological diversity we carried out a pixel-based analysis using a linear modelling framework (Faraway 2014).In the first instance, we used linear modelling framework to quantify changes in phenological diversity between the two periods, 2003-2010 and 2011-2020.The model takes the form: where: X i 1 b ( ) is a two-level categorical variable representing the two periods (2003-2010 and 2011-2020)   and: N 0, 2 e s ~( ) We further quantified sensitivity of changes in phenology diversity to two key climatic drivers, temperature and cumulated precipitation.Time series data were derived from the Terraclimate dataset (Abatzoglou et al 2018).
The model takes the form: where: X 1 b is a vector temperature or precipitation data and: ) We filtered out non-statistically significant results and retained those where the P value for the parameter estimates (P < 0.05).

Drivers of spatial patterns in phenological diversity
Spatial patterns in phenological diversity show significant differences among biomes, with temperate and arid regions displaying the highest levels of phenological diversity (figure 2).On the contrary, the boreal and tropical biomes show the lowest levels of phenological diversity.A variance component analysis of the spatial patterns of variability (see Methods) among biomes and plant functional types, revealed that among-biome variability contributed to 18.76% of the total variation in phenological diversity.The component of variation uniquely attributable to plant functional types is small, at only 1.83%.The largest variance component is the residual error (or the population-level variance), at 79.41% of the total variance.The Random Forests analysis revealed that climatic variables had a greater influence than topographic or the human footprint variables.Specifically, temperature seasonality and mean annual temperature were of greatest importance in explaining phenological diversity (figure 3(A)).The influence of other variables ranged from moderate to low (figure 3(A)).These included spatial heterogeneity in mean annual temperature and aridity, which were moderately influential in the model (ranked 4th and 5th, respectively, see figure 3(A)).However, their importance was only marginally higher compared to other variables, i.e. precipitation seasonality, standard deviation (hereafter sd) in temperature seasonality, and annual precipitation (sd Temperature Seasonality and sd Precipitation, see figure 3(A)). .Overall, the ten predictors explained 43% of the spatial variability in phenological diversity, as obtained from the cross-validation procedure (figure S2).
An inspection of the functional form of the fitted relationships, between phenological diversity and individual predictors (figure S3), revealed that responses were generally non-linear.These included a broadly positive effect of variables related to spatial heterogeneity in climate (sd Annual Precipitation, sd Temperature seasonality, sd Annual mean temperature).Other variables, such as aridity and mean annual temperature, showed a more complex response.This complexity was further confirmed by a detailed analysis of the local effects of the predictors, which highlighted marked differences in responses across biomes (figure S4 and figure S5).Notably, in the boreal biome, responses to mean annual temperature and temperature seasonality were primarily negative.On the contrary, in other biomes, these variables had an overwhelmingly positive effect.Among the water-related variables, aridity had a predominantly negative effect across the boreal, temperate, and arid biomes.
Bivariate partial dependence plots, depicting the effects of interactions between predictors, were also important in assessing the drivers of phenological diversity (figures 3(B)-(F)).Complex interactions were .Specifically, a combination of annual mean temperatures above ∼0 °C and high levels of aridity, had a generally positive effect on phenological diversity (figure 3(C)).Conversely, mean annual temperatures below ∼0 °C had a negative impact on phenological diversity, regardless of aridity levels.Seasonality variables showed a positive interactive effect when temperature seasonality was below ∼ 8 °C.On the contrary, temperature seasonality values above this threshold showed a negative interaction with precipitation seasonality (figure 3(C)).Less complex interactions were detected for other variables, describing spatial heterogeneity in climate, topography and the human footprint index, suggesting that additive or linear effects are likely more important for these variables (figures 3(D)-(F)).

Temporal changes in phenological diversity
The availability of a 20-year long data series from the MODIS sensor, combined with the high temporal consistency of retrievals derived from this platform, supports the assessment of the emerging recent trends in phenological diversity.Linear regression models revealed that phenological diversity had changed (P < 0.05) over 21% of the analysed area, when comparing the two periods 2003-2011 and 2012-2020.A summary of the model coefficients indicates that, on average, the strongest reduction in phenological diversity occurred within the boreal biome (figure 4(A)).The temperate biome also showed a reduction in phenological diversity, albeit more modest.On the contrary, arid and tropical biomes showed, on average, an increase in phenological diversity (figure 4(A)).A summary of the statistically significant results of the linear models (P < 0.05), aimed at quantifying the sensitivity of phenological diversity to temperature and precipitation, revealed that 12% and 11% of the coefficients were statistically significant, respectively.In general, we diagnose that warming has led primarily to a negative effect on phenological diversity within the boreal and temperate biomes (figure 4(B)).However, the opposite was true for the arid and tropical biomes, where increasing temperature positively affected phenological diversity (figure 4(B)).Models including precipitation as a predictor showed that, on average, precipitation had a negative effect on phenological diversity in temperate and arid biomes (figure 4(C)).In contrast, the effect of precipitation had mostly a positive effect in the boreal and tropical biomes (figure 4(C)).

Discussion
Phenology is a key plant functional trait that can be effectively retrieved from Earth Observation data.Quantifying spatial patterns of heterogeneity in phenological strategies can thus offer significant promise for advancing under-explored areas of plant functional ecology, beyond vegetation-climate relationships alone (Viña et al 2012, 2016, Dronova and Taddeo 2022).
Phenological diversity is not only thought to reflect the functional composition of ecosystems, but also their resilience, because of its ability to inform on the level of synchronicity in the ecosystem demand for environmental resources (Lasky et al 2016, Valencia et al 2020).
Our results showed important differences in phenological diversity patterns among biomes, which may depend on the (co-)occurrence of different phenological strategies (Harrison et al 2010).The low levels of phenological diversity observed in boreal and humid tropical forests is likely due to the abundance of evergreen tree species with a homogeneous phenological strategy, including evergreen broadleaf species in humid tropical forests, and evergreen needle-leaf species in boreal areas.On the contrary, temperate and arid biomes, which show the highest levels of phenological diversity, host a number of co-occurring tree species, with different phenological strategies.These include species in mixed deciduous /evergreen and deciduous stands in temperate areas, and evergreen and deciduous species in arid areas.The high levels of phenological diversity within these biomes could be due to a number of non-mutually exclusive factors.For instance, in sparse forests the signal could be affected by understory species with different phenological strategies, which are likely to have an impact on the emergent phenological diversity.
In arid biomes, the presence of species with alternative photosynthetic pathways, combined with different strategies for the use of water resources, may result in a high level of phenological diversity (Biederman et al 2017, Smith et al 2019, Kato et al 2021).It is also possible that, in some cases, mixed pixels with a gradient of forest to grass cover could be influencing the observed patterns of phenological diversity.Ultimately, the presence of hotspots of phenological diversity within the temperate and arid biomes, could be due to a number of ecophysiological constraints, including a co-limitation of primary productivity, driven by different factors including water availability and temperature (Harrison et al 2010, Huang et al 2017, Peaucelle et al 2019).These factors can lead to higher spatial variability in environmental conditions, which in turn translates into higher phenological diversity.
Our Random Forests analysis showed that there was a strong climatic control on phenological diversity.The role of climate as a driver of phenology is well documented (Richardson et al 2013, Forkel et al 2015), but this had not been tested for the spatial heterogeneity of this functional property.This finding further reinforces the usefulness of phenology as an indicator of forest ecosystem responses to environmental change.
Temperature-related variables appear to be primary predictors of phenological diversity.However, three of the water-related variables contributed to over 29% of the explained variance, suggesting a joint role of temperature and precipitation in shaping phenological diversity patterns.
We note that the phenological diversity of boreal and temperate regions tends to be more related to temperature-related variables (Supplementary figure 4 and figure 5).On the contrary, in arid and tropical regions, phenological diversity is more controlled by water-related variables.These nonstationary, i.e. spatiallyvarying, effects suggest that different factors may drive phenological diversity at regional scales.We interpret these different responses to climatic drivers at regional scales as a reflection of the response to the local limiting factors of the dominant tree species, which in turn may exert a strong influence on ecosystem functional properties, via the mass-ratio effect (Grime 1998).
The analysis of the annual time series data provided two key insights into spatio-temporal patterns of change in phenological diversity.First, our results show that the boreal, temperate and arid biomes are areas of rapid ecological change, with marked variations in phenological diversity.Interestingly, biomes in the northern and southern hemispheres show a contrasting direction of change, with the boreal and temperate biomes displaying primarily a decline in phenological diversity and the tropical and arid biomes showing an increase in phenological diversity.These findings suggest that a certain degree of reshuffling in phenological strategies might be occurring, as part of a generalised trend in the changes in the seasonal profile of plant vegetation growth across these biomes (Buitenwerf et al 2015, Forkel et al 2015, Garonna et al 2016, Huang et al 2017).
Second, changes in phenological diversity exhibit biome-specific sensitivity to fluctuations in temperature and precipitation regimes.We tentatively attribute these differences in sensitivities to factors controlling seasonal plant growth.For instance, warming in the boreal zone may lead to an enhancement/anticipation of phenology.On the contrary, arid regions show the opposite effect with a decrease in phenological diversity, when warming occurs.These sensitivities to temperature and precipitation might be further modulated by differences in the distribution of tree phenological strategies (e.g.evergreen versus deciduous) across biomes.
Monitoring changes in the functional composition of forest ecosystems requires robust metrics that can be compared across large areas.The metric developed here has the advantage of minimizing the temporal white noise of satellite retrievals by deriving phenological information from the cumulated sum of a spectral index.Furthermore, our metric is based on the variation of phenological trajectories around the average phenology for a given area, thus making it insensitive to mean local phenological patterns, and enabling a comparison of phenological diversity across biomes.
Importantly, our metric can be applied at a global scale using data from different satellite platforms.We envisage that the use of data with a higher spatial resolution, such as those generated by the ESA Sentinel 2 satellites (Berger et al 2012, Malenovský et al 2012) of the Copernicus programme, will significantly improve our understanding of the relationship between phenological diversity and the functional structure of tree communities through the fine-scale mapping capabilities of novel sensors.
While land-surface phenology has primarily been advocated as a metric for monitoring ecosystem responses to environmental change, we show that it has also the potential to be used as an indicator of the functional composition of ecosystems.
As a trait that scales up across different levels of biological organisation, phenology can also be interpreted in the light of processes operating at specific scales.As recently discussed by Dronova and Taddeo (2022), phenological diversity metrics computed from land surface phenology have the potential to provide key insights into ecological processes occurring at the level of tree communities.Indicators of phenological diversity may thus be ideal candidates for Essential Biodiversity Variables, potentially adding a new group of metrics to the current portfolio of proposed ecosystem-level Essential Biodiversity Variables (Skidmore et al 2021).
There is an increasing need for monitoring plant functional traits to assess land-climate interactions and the impacts of climate change on ecosystems (Violle et al 2014).To this end, global information at high spatial resolution is required to drive climate adaptation policies (Commission 2020b, 2020a).Our work ultimately contributes to the development an operational framework for quantifying spatio-temporal changes in ecosystem functional composition by emphasising the complex interactions between background climate conditions, climate change and phenological diversity of the terrestrial biosphere.

Figure 1 .
Figure 1.Phenological diversity metric algorithm workflow (see methods for a detailed description of each step).

Figure 2 .
Figure 2. Global patterns of forest phenological diversity.The histogram on the bottom-left corner shows the proportion of the pixels that fall within each of the discrete colour classes.The left-panel shows biome-level summaries in the variation of phenological diversity.Percentiles (25th, 50th, 75th) are shown as vertical bars within the empirical distributions.The 50th percentile i.e. the median is shown in red.

Figure 3 .
Figure 3. Relative importance of the predictor variables for the random forest regression model (panel A).Colours distinguish the different categories of environmental predictors (climate, topography and human impact).Two-dimensional partial dependence plots showing the interactions between predictors and their effects on phenological diversity (panels B-F).All variables except those graphed are held at their means.

Figure 4 .
Figure 4. Biome-level changes in phenological diversity between two periods 2011-2020 and 2003-2010 (A), sensitivity to temperature (B) and precipitation (C).Changes and sensitivity in phenological diversity are quantified by summarizing the pixel-level effects estimated via linear modelling.Coloured circles indicate the median response and triangles show the 25th and 75th percentile for each biome.Distance from zero of the distribution is showed by an horizontal bar.The red line shows the global average.