Phenological Characteristics of Global Ecosystems Based on Optical, Fluorescence, and Microwave Remote Sensing

: Growing seasons of vegetation generally start earlier and last longer due to anthropogenic warming. To facilitate the detection and monitoring of these phenological changes, we developed a discrete, hierarchical set of global “phenoregions” using self-organizing maps and three satellite-based vegetation indices representing multiple aspects of vegetation structure and function, including the normalized di ﬀ erence vegetation index (NDVI), solar-induced chlorophyll ﬂuorescence (SIF), and vegetation optical depth (VOD). Here, we describe the distribution and phenological characteristics of these phenoregions, including their mean temperature and precipitation, di ﬀ erences among the three satellite indices, the number of annual growth cycles within each phenoregion and index, and recent changes in the land area of each phenoregion. We found that the phenoregions “self-organized” along two primary dimensions: degree of seasonality and peak productivity. The three satellite-based indices each appeared to provide unique information on land surface phenology, with SIF and VOD improving the ability to detect distinct annual and subannual growth cycles in some regions. Over the nine-year study period (limited in length by the short satellite SIF record), there was generally a decrease in the spatial extent of the highest productivity phenoregions, though whether due to climate or land use change remains unclear.


Introduction
The timing of vegetation phenology-the seasonal cycles of vegetation growth, development, and reproduction-depends on cues from environmental conditions like daylength, temperature, and moisture availability [1][2][3]. Warming has generally led to an earlier start and an overall lengthening of the vegetation growing season, observed both in situ and with remote sensing of "land surface phenology" [2][3][4], the seasonal variation of land surface "greenness" observed at synoptic scale [5]. Phenological change can also result from nonclimatic forcings, like land use and land cover change [6][7][8] and disturbance [9], which should be controlled when attributing observed phenology changes to climate. Detecting and monitoring phenological responses to climate variability and climate change with remote sensing is simplified when performed on "phenoregions" [7,10] (clusters of pixels with self-similar phenological and climatic characteristics), especially when excluding phenoregions with significant nonclimatic changes [7]. Rather than analyzing individual pixels, phenoregion-scale monitoring improves signal retrieval by reducing the influence of noise, facilitating interpretation of change, and eliminating the need for forward-or backward-looking smoothing functions [10]. The delineation of phenoregions at regional and global scales is typically based on clustering time series of satellite-derived normalized difference vegetation index (NDVI) along with other ancillary datasets (e.g., gridded climate data) [7,11,12].
Although the normalized difference vegetation index (NDVI) may perform well as the sole phenological indicator in temperate, deciduous-dominated ecosystems, it is not well suited for phenoregion delineation in regions with prolonged durations of cloud cover (e.g., the tropics), with predominantly evergreen leaf phenologies, or with high spatiotemporal heterogeneity of land surface phenology (e.g., in dryland ecosystems) [7,13]. Since NDVI primarily responds to canopy greenness dynamics resulting from changes in leaf area and chlorophyll content [14], it also may not capture the full range of functional and structural manifestations of phenology (e.g., leaf area phenology, leaf pigment phenology, photosynthetic phenology, reproductive phenology, etc.), especially where physiological function is decoupled from canopy greenness (e.g., in the tropics and in drylands) [3,13]. For these reasons, the tropics and many dryland regions are often excluded from phenoregion classifications due to lack of an observable annual NDVI cycle [7]. To improve understanding and monitoring of phenological dynamics outside the temperate and boreal regions, alternative vegetation indices are necessary to capture a greater range of ecosystem structural and functional characteristics.
Two emerging remotely sensed indicators could be particularly useful for capturing structural and functional dynamics that are not available from NDVI alone: microwave-based vegetation optical depth (VOD) and solar-induced chlorophyll fluorescence (SIF). VOD, measuring microwave signal attenuation through plant canopies, is capable of tracking vegetation water content in both photosynthetic and non-photosynthetic tissues of aboveground biomass [15] while being less sensitive to cloud and aerosol contamination [16,17], making it complementary to NDVI for tracking changes in woody biomass [18] and in cloud-prone regions [15,17,19,20]. Whereas NDVI and VOD both track relatively slow-varying phenological processes (i.e., development of the canopy and accumulation of biomass), chlorophyll fluorescence (the reemission of photons by chlorophyll during photosynthesis) can capture the fast-changing dynamics of plant physiological function [21][22][23]. SIF can therefore better track vegetation dynamics in regions where canopy greenness is decoupled from plant physiological activity, particularly in evergreen forests [21,24] and dryland ecosystems [25,26].
Here, we delineated global, coarse-resolution phenoregions based on NDVI, SIF, and VOD using a self-organizing map (SOM) method [27]. Using these phenoregions, we explored multiple characteristics of global vegetation phenology, including: (1) Differences in phenology among the global phenoregions and among the remotely sensed indicators; (2) the climatic characteristics of the phenoregions; (3) the presence or absence of (sub)annual growing season dynamics within each phenoregion and vegetation index; and (4) recent trends in the global land area of each phenoregion. As these phenoregions have relatively coarse resolutions (0.5 • grid cells) due to limitations of the SIF and VOD data, this research highlights the utility of the SOM method for phenoregion delineation and the promise of using a broad suite of remotely sensed data (complementing traditional greenness indices) to capture multiple structural and functional aspects of land surface phenology.

Remotely Sensed Vegetation Information
We obtained NDVI, SIF, and VOD data over a common period of overlap (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). The datasets include twice-monthly, 8-km resolution NDVI from the Global Inventory Modeling and Mapping Studies (GIMMS) [28]; daily 25-km resolution X-band (10.7 GHz) VOD from the Advanced Microwave Scanning Radiometer (AMSR-E) on NASA's EOS-Aqua satellite and the Advanced Microwave Scanning Radiometer 2 (AMSR2) on the JAXA GCOM-W1 satellite [29]; and monthly 0.5 • resolution Version 2.6 SIF from the Global Ozone Mapping Experiment (GOME-2) [30]. The NDVI and VOD data were averaged to monthly, 0.5 • resolution to match the SIF data. Before resampling, poor quality pixels were Remote Sens. 2020, 12, 671 3 of 15 filtered from the three datasets. NDVI values were removed when flagged as snow or cloud covered. We also excluded VOD values affected by frozen ground, snow or ice, strong precipitation, or large retrieval uncertainty based on the quality flag provided by the Version 2 global land parameter data record [15,29]. The SIF data were filtered to exclude observations with cloud cover or with solar zenith angle >70 • [30,31]. To ensure that all datasets were on a similar scale prior to clustering, we divided the VOD (which is in arbitrary units) by two.

Defining Hierarchical Phenoregions Using Self-Organizing Maps
To characterize global land surface phenology, we grouped pixels into phenoregions based on the mean seasonal cycles of NDVI, SIF, and VOD. The clusters were defined using self-organzing maps (SOM), a neural network method where neurons (or "nodes")-each associated with an n-dimensional "model" vector (where n is the number of dimensions, or variables, in the dataset)-are arranged in a grid, or "map", representing the data space in fewer dimensions, typically a one or two dimensional array ( Figure 1) [27]. The nodes are topologically organized such that nodes that are closer together in the SOM space are also more similar in the n-dimensional hyperspace defined by the data, with adjacent nodes connected by a neighborhood function. The nodes in the SOM can be thought of as clusters to which each observation in the dataset can be assigned based on the minimum Euclidean distance between the n-dimensional observations and the model vectors of each node.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 14 as snow or cloud covered. We also excluded VOD values affected by frozen ground, snow or ice, strong precipitation, or large retrieval uncertainty based on the quality flag provided by the Version 2 global land parameter data record [15,29]. The SIF data were filtered to exclude observations with cloud cover or with solar zenith angle >70° [30,31]. To ensure that all datasets were on a similar scale prior to clustering, we divided the VOD (which is in arbitrary units) by two.

Defining Hierarchical Phenoregions Using Self-Organizing Maps
To characterize global land surface phenology, we grouped pixels into phenoregions based on the mean seasonal cycles of NDVI, SIF, and VOD. The clusters were defined using self-organzing maps (SOM), a neural network method where neurons (or "nodes")-each associated with an ndimensional "model" vector (where n is the number of dimensions, or variables, in the dataset)-are arranged in a grid, or "map", representing the data space in fewer dimensions, typically a one or two dimensional array (Figure 1) [27]. The nodes are topologically organized such that nodes that are closer together in the SOM space are also more similar in the n-dimensional hyperspace defined by the data, with adjacent nodes connected by a neighborhood function. The nodes in the SOM can be thought of as clusters to which each observation in the dataset can be assigned based on the minimum Euclidean distance between the n-dimensional observations and the model vectors of each node. Whereas, in some respects, the SOM is similar conceptually to other clustering methods (e.g., kmeans clustering), it has the benefit of a clear organizational structure of the clusters, thus allowing discretization of input data (e.g., a discrete set of phenoregions) while explicitly preserving a degree of continuity among the clusters. This property has led to the increasingly widespread use of the SOM as a clustering method in the earth sciences, particularly in synoptic climatology, and, more recently, in the remote sensing of vegetation. Within climatology, the SOM is often used to characterize dominant spatial patterns of climate (in which each "node" in the SOM represents a synoptic-scale spatial pattern) that vary through time [32,33]. Self-organizing maps have thus been used to delineate dominant patterns of sea level pressure [34], 500 mb geopotential heights [35], sea surface temperature [36], and surface weather patterns [37]. SOMs have also been used within remote sensing to characterize dominant spatial patterns of vegetation index anomalies [38]. In these cases, the observational units being clustered are spatial patterns at particular points in time, whereas in the case of a phenological classification, the observations and variables are reversed, with pixels as the units being classified on the basis of their temporal patterns. Whereas, in some respects, the SOM is similar conceptually to other clustering methods (e.g., k-means clustering), it has the benefit of a clear organizational structure of the clusters, thus allowing discretization of input data (e.g., a discrete set of phenoregions) while explicitly preserving a degree of continuity among the clusters. This property has led to the increasingly widespread use of the SOM as a clustering method in the earth sciences, particularly in synoptic climatology, and, more recently, in the remote sensing of vegetation. Within climatology, the SOM is often used to characterize dominant spatial patterns of climate (in which each "node" in the SOM represents a synoptic-scale spatial pattern) that vary through time [32,33]. Self-organizing maps have thus been used to delineate dominant patterns of sea level pressure [34], 500 mb geopotential heights [35], sea surface temperature [36], and surface weather patterns [37]. SOMs have also been used within remote sensing to characterize dominant spatial patterns of vegetation index anomalies [38]. In these cases, the observational units being clustered are spatial patterns at particular points in time, whereas in the case of a phenological classification, the observations and variables are reversed, with pixels as the units being classified on the basis of their temporal patterns. Defining the number of clusters (or "nodes") to use, as well as their arrangement in the low-dimensional SOM space, is the key decision that needs to be made [27,32,33]. Although there are some objective criteria available to make this decision [e.g., 35,36], there is no clear consensus on the best approach, and many clustering applications set relatively arbitrary numbers of clusters based on judgement of the investigator, ultimately attempting to balance competing needs for accuracy, specificity, and ease of visualization [27]. Using the SOM Toolbox in MATLAB [39], we chose to emphasize ease of visualization in a simple proof-of-concept analysis, with a small number of Level I clusters (12 total) in a rectangular three-by-four SOM arrangement ( Figure 1) based on 36 variables: mean monthly NDVI, SIF, and VOD over the period 2007-2015. Maintaining this small number of Level I phenoregions has the added benefit of being roughly consistent with the number of regions in the Köppen climate classification system and the EPA Level I ecoregion classification, which could allow more direct comparison with these other products. To maintain consistent growing season timing across the northern and southern hemispheres, we used mean monthly vegetation indices from January to December for the Northern Hemisphere and from July to June for the Southern Hemisphere.
After defining the Level I clusters, we used additional SOM models to further subdivide each cluster into six additional Level II clusters based on the same 36 remotely sensed variables, thus revealing additional phenological detail within the Level I clusters. Since the SOMs are arranged in two-dimensional arrays, the Level I SOMs are constrained to only represent the two primary dimensions of the data, so subdividing each phenoregion by a second SOM will effectively allow "rotation" onto other dimensions of the data (e.g., number of annual growing seasons or timing of growing seasons). We chose to divide each Level I phenoregion into six Level II regions (in a two-by-three rectangular arrangement) to maintain class sizes well in excess of Kohonen's recommendation of at least 50 observations per SOM node [27].

Phenological Characteristics of the Global Phenoregions
We examined the climatic characteristics of each phenoregion using 0.5 • monthly mean temperature and precipitation from the University of East Anglia's Climatic Research Unit (CRU) TS4.01 product [40]. Specifically, we extracted temperature and precipitation estimates for each phenoregion and calculated the mean (and standard deviation) of temperature and precipitation for each month based on the period of overlap with the satellite data (2007-2015).
We tested for annual (one cycle per year) and biannual (two cycles per year) phenological cycles in each phenoregion and for each vegetation index using the multitaper method (MTM), a robust approach for estimating power spectral density (PSD) [41]. Annual and biannual signals are not mutually exclusive and any given pixel may contain elements of both [42], often with a lower-amplitude biannual signal overlain on a dominant annual signal (e.g., Figure 2). At each 0.5 • pixel, we calculated the PSD on the full time series of each vegetation index (12 months × 9 years = 108 observations per index per pixel) and examined the PSD for significant signals at annual and biannual frequencies. The significance of the PSD at a given frequency was assessed by comparing the observed PSD to 1000 realizations of pink noise with the same mean and variance as the observations (Figure 2), testing the null hypothesis that the observed PSD did not exceed that expected from random noise, with the null hypothesis rejected if the observed PSD exceeded the 99 th percentile of the PSD from the noise-only ensemble. Since the MTM requires a complete time series (i.e., no missing observations), we filled missing values first by replacing them with the pixel's seasonal mean values. Remaining missing values (i.e., where monthly means could not be defined due to persistent missing values) were filled using linear interpolation, which was particularly common in the VOD data during winter months at high latitudes and elevations. Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 14 Finally, we examined recent trends in the land area of each phenoregion by annually reassigning each pixel to its closest matching phenoregion based on the annual NDVI, SIF, and VOD observations. This allowed pixels to move in and out of phenoregions based on interannual changes in the land surface phenology observed by the satellites. We first calculated the nearest Euclidean distance between the model vector of each SOM node (i.e., phenoregion) and the monthly NDVI, SIF, and VOD observations for each individual year and pixel. Here, we used only quality-screened vegetation index observations, with missing or low-quality observations excluded from the distance calculation to avoid introducing uncertainties from the gap filling process. After annually reassigning each pixel to its closest-matching phenoregion, we calculated the total annual land area (in km 2 ) assigned to each phenoregion (correcting for the effect of latitude on 0.5° pixel size) and assessed trends in the phenoregion land area during the study period (2007-2015) using ordinary least squares regression.

Distribution and Characteristics of Global Phenoregions
Global ecosystems "self-organized" into clusters differentiated predominantly along two main phenological dimensions (Figure 3): degree of seasonality (ranging from high seasonality in the first row of the SOM to low seasonality in the bottom row) and peak productivity (ranging from low productivity in the left-hand columns of the SOM to high productivity in the right-hand column). The most seasonal phenoregions (phenoregions 1-4) were largely confined to the mid-to-high latitudes of the Northern Hemisphere ( Figure 4), with very large annual temperature ranges and low precipitation (well below 100 mm/month in most regions) ( Figure 5). The phenoregions with the Finally, we examined recent trends in the land area of each phenoregion by annually re-assigning each pixel to its closest matching phenoregion based on the annual NDVI, SIF, and VOD observations. This allowed pixels to move in and out of phenoregions based on interannual changes in the land surface phenology observed by the satellites. We first calculated the nearest Euclidean distance between the model vector of each SOM node (i.e., phenoregion) and the monthly NDVI, SIF, and VOD observations for each individual year and pixel. Here, we used only quality-screened vegetation index observations, with missing or low-quality observations excluded from the distance calculation to avoid introducing uncertainties from the gap filling process. After annually reassigning each pixel to its closest-matching phenoregion, we calculated the total annual land area (in km 2 ) assigned to each phenoregion (correcting for the effect of latitude on 0.5 • pixel size) and assessed trends in the phenoregion land area during the study period (2007-2015) using ordinary least squares regression.

Distribution and Characteristics of Global Phenoregions
Global ecosystems "self-organized" into clusters differentiated predominantly along two main phenological dimensions (Figure 3): degree of seasonality (ranging from high seasonality in the first row of the SOM to low seasonality in the bottom row) and peak productivity (ranging from low productivity in the left-hand columns of the SOM to high productivity in the right-hand column). The most seasonal phenoregions (phenoregions 1-4) were largely confined to the mid-to-high latitudes of the Northern Hemisphere (Figure 4), with very large annual temperature ranges and low precipitation (well below 100 mm/month in most regions) ( Figure 5). The phenoregions with the highest seasonality and productivity (phenoregions 4 and 8) included much of North America's and Eurasia's croplands (Figures 3 and 4), which were largely distinguished by their high growing season SIF [43].
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 14 highest seasonality and productivity (phenoregions 4 and 8) included much of North America's and Eurasia's croplands (Figures 3 and 4), which were largely distinguished by their high growing season SIF [43].  The legend is organized in the same two-dimensional arrangement as the SOM (with mean seasonal cycles shown in Figure 3), with increasing peak productivity from left to right (shaded from red to purple) and increasing seasonality from bottom to top (shaded from darker to lighter). To maintain consistent growing season timing across hemispheres, the seasonal cycles shown above start in January for the Northern Hemisphere and July for the Southern Hemisphere.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 14 highest seasonality and productivity (phenoregions 4 and 8) included much of North America's and Eurasia's croplands (Figures 3 and 4), which were largely distinguished by their high growing season SIF [43].  The legend is organized in the same two-dimensional arrangement as the SOM (with mean seasonal cycles shown in Figure 3), with increasing peak productivity from left to right (shaded from red to purple) and increasing seasonality from bottom to top (shaded from darker to lighter). The legend is organized in the same two-dimensional arrangement as the SOM (with mean seasonal cycles shown in Figure 3), with increasing peak productivity from left to right (shaded from red to purple) and increasing seasonality from bottom to top (shaded from darker to lighter).  The least seasonal phenoregions (phenoregions 9-12; Figure 3) were mostly distributed across global tropical and arid regions (Figure 4). The least productive of these regions (phenoregions 9-10) tended to have very low precipitation (<100 mm/month) and highly variable seasonal temperatures ( Figure 5). The most extensive phenoregion (phenoregion 9, covering 24% of global land area; Figure  4) largely consisted of relatively unvegetated hyperarid deserts, including the Sahara, Arabian, Namib, and Atacama deserts, as well as much of central Asia and Australia. The more productive phenoregions (phenoregions 11-12) tended to have high temperatures throughout the year and high precipitation, especially in summer ( Figure 5). These Level I phenoregions clearly captured rainfall gradients, most clearly in the transition from the Sahara Desert (phenoregion 9) to the Sahel (phenoregions 10-11) to the Congo rainforests (phenoregion 12).
Even in the least seasonal phenoregions, the mean SIF seasonal cycle showed substantial intraannual phenology compared to the mean NDVI and VOD seasonal cycles, especially in the tropics and other high-productivity regions (Figure 3). Reflectance-based indices like NDVI largely track morphological development of the canopy (e.g., green leaf area) rather than physiological processes, so NDVI remained relatively constant in these evergreen systems, but SIF appeared to track their photosynthetic phenology (Figure 3) [44]. Likewise, in phenoregions 10 and 11, which generally occurred in semiarid and subtropical regions, the mean SIF seasonal cycle appeared to track both spring and late summer growing season peaks that did not appear in the NDVI or VOD data. This provides additional evidence that SIF may be capable of tracking intra-annual photosynthetic phenological dynamics that are not tracked by traditional multispectral vegetation indices [25,45], especially in semiarid systems [13,25,[46][47][48][49][50][51][52].
Whereas the Level I phenoregions were predominantly separated based on seasonality and peak productivity, the Level II phenoregions (Supplemental Figures S1-S12) were further distinguished by growing season timing and number of growing seasons per year. Within phenoregion 11, for example, which largely occurred along the periphery of (sub)tropical forests and transition zones between the tropics and semiarid regions, the Level II phenoregions clearly separated regions dominated by summer and winter rains ( Figure S11). The Sahel dominated a cluster (phenoregion 11.4) consisting of a single-peaked growing season centered in late summer. This phenoregion was differentiated from others with winter/spring growing seasons (characteristic of Mediterranean The least seasonal phenoregions (phenoregions 9-12; Figure 3) were mostly distributed across global tropical and arid regions (Figure 4). The least productive of these regions (phenoregions 9-10) tended to have very low precipitation (<100 mm/month) and highly variable seasonal temperatures ( Figure 5). The most extensive phenoregion (phenoregion 9, covering 24% of global land area; Figure 4) largely consisted of relatively unvegetated hyperarid deserts, including the Sahara, Arabian, Namib, and Atacama deserts, as well as much of central Asia and Australia. The more productive phenoregions (phenoregions 11-12) tended to have high temperatures throughout the year and high precipitation, especially in summer ( Figure 5). These Level I phenoregions clearly captured rainfall gradients, most clearly in the transition from the Sahara Desert (phenoregion 9) to the Sahel (phenoregions 10-11) to the Congo rainforests (phenoregion 12).
Even in the least seasonal phenoregions, the mean SIF seasonal cycle showed substantial intra-annual phenology compared to the mean NDVI and VOD seasonal cycles, especially in the tropics and other high-productivity regions (Figure 3). Reflectance-based indices like NDVI largely track morphological development of the canopy (e.g., green leaf area) rather than physiological processes, so NDVI remained relatively constant in these evergreen systems, but SIF appeared to track their photosynthetic phenology (Figure 3) [44]. Likewise, in phenoregions 10 and 11, which generally occurred in semiarid and subtropical regions, the mean SIF seasonal cycle appeared to track both spring and late summer growing season peaks that did not appear in the NDVI or VOD data. This provides additional evidence that SIF may be capable of tracking intra-annual photosynthetic phenological dynamics that are not tracked by traditional multispectral vegetation indices [25,45], especially in semiarid systems [13,25,[46][47][48][49][50][51][52].
Whereas the Level I phenoregions were predominantly separated based on seasonality and peak productivity, the Level II phenoregions (Supplemental Figures S1-S12) were further distinguished by growing season timing and number of growing seasons per year. Within phenoregion 11, for example, which largely occurred along the periphery of (sub)tropical forests and transition zones between the tropics and semiarid regions, the Level II phenoregions clearly separated regions dominated by summer and winter rains ( Figure S11). The Sahel dominated a cluster (phenoregion 11.4) consisting of a single-peaked growing season centered in late summer. This phenoregion was differentiated from others with winter/spring growing seasons (characteristic of Mediterranean ecosystems in California and Spain; phenoregion 11.3) and those with a primary summer growing season but a secondary peak in winter/spring (phenoregion 11.1). Likewise, phenoregion 12, largely consisting of tropical and subtropical forests, was further split by the number, timing, and magnitudes of annual growth cycles ( Figure S12). The equatorial tropics largely fell into one of two Level II phenoregions: One with two roughly equal SIF peaks in spring and late summer, but very low seasonality in NDVI and VOD (phenoregion 12.3), and one with a biannual SIF cycle overlain on a distinct annual cycle with small winter peaks in VOD and NDVI (phenoregion 12.6). These were distinct from phenoregion 12.4, which featured a single-peaked SIF phenological cycle (with maxima in summer) and smaller annual cycles in NDVI and VOD, largely occurring around the margins of equatorial tropical forests and in the southeastern United States, western Europe, and southeast Asia.
This initial work demonstrates both the utility of the SOM method and the promise of incorporating multiple vegetation indices, but future work could further refine and improve these phenoregions. First, while SIF shows promise for capturing photosynthetic phenologies in many systems, it is limited by a short temporal record (only back to 2007 for GOME-2) and coarse spatial resolution (a 0.5 • grid for GOME-2) [3]. With the advent of new and improved SIF sensors, such as OCO-3 [53] and TROPOMI [54], it may soon be possible to improve the spatial resolution of SIF-based phenoregions, which would be especially useful for applications focused on separating climate from land use drivers of phenological change, since land use change typically occurs at a much finer scale than the 0.5 • resolution used here. Second, we chose the number of phenoregions primarily to prioritize ease of visualization. Using more objective methods (e.g., [35,36]) to formally define the number of statistically distinguishable phenoregions could provide a more complete and well-constrained phenoregion map.

Spectral Analysis of Phenological Dynamics and Differences among Vegetation Indices
Globally, most pixels had distinct annual cycles in all three vegetation indices (Figures 6 and 7), but the regions and indices differed regarding the presence of biannual seasonal dynamics and in the spatial patterns of these (bi)annual cycles. In the most seasonal phenoregions (phenoregions 1-4), the MTM power spectral density detected significant annual phenologies in 90%-100% of pixels, but SIF only detected an annual cycle in~80% of pixels in the lowest productivity systems (phenoregion 1). However, the SIF MTM power spectrum also showed significant biannual phenological signals overlain on the annual cycles (e.g., Figure 2) in~20% of pixels (largely in northern and western Eurasia; Figure 7) that were far less prevalent in the NDVI and VOD, suggesting that SIF may detect subannual photosynthetic phenologies that are not present in canopy morphology or water content of aboveground biomass.
The MTM power spectrum detected fewer distinct growing seasons for all three remotely sensed indices in tropical and arid regions (Figures 6 and 7). In the Amazon, VOD generally detected more distinct annual growing seasons than NDVI or SIF ( Figure 7A-C), possibly reflecting the relatively low sensitivity of VOD to the persistent cloud cover in this region [15,17,19,20]. All three indices, especially SIF and NDVI, captured clear biannual phenological dynamics in tropical forests of Africa ( Figure 7D-F). SIF and VOD also showed both annual and biannual growing seasons in parts of the Sahel, whereas NDVI did not clearly detect the biannual component of the growing season in this region ( Figure 7D-F).
In the least productive regions (phenoregions 5 and 9), SIF detected fewer distinct growing seasons than NDVI or VOD (Figures 6 and 7). This may be somewhat surprising given the presence of clearer seasonal cycles in the mean SIF phenology (Figure 3), but we suspect that this discrepancy can be traced to one of two possible explanations. First, averaging across space and time (as in Figure 3) likely cancelled out some of the relatively high random noise in the SIF signal, leaving a more visually obvious phenological cycle than that achieved from the full monthly time series used in the MTM, where the satellite-based SIF signal would be less distinct from the noise spectrum. Alternatively, in the lowest productivity phenoregions (mostly arid or hyperarid regions), vegetation cover is exceptionally low, and any phenological cycle observed at the satellite is likely unrelated to phenological changes in the morphology or physiology of vegetation. The lack of obvious seasonal cycle in the SIF power spectrum may therefore more closely approximate the true vegetation cycles in these regions, while the apparent cycles seen in the NDVI and VOD could be artifacts of sun-sensor geometry, retrieval error, variations in soil color/wetness, or atmospheric effects.

Spectral Analysis of Phenological Dynamics and Differences among Vegetation Indices
Globally, most pixels had distinct annual cycles in all three vegetation indices (Figures 6 and 7), but the regions and indices differed regarding the presence of biannual seasonal dynamics and in the spatial patterns of these (bi)annual cycles. In the most seasonal phenoregions (phenoregions 1-4), the MTM power spectral density detected significant annual phenologies in 90%-100% of pixels, but SIF only detected an annual cycle in ~80% of pixels in the lowest productivity systems (phenoregion 1). However, the SIF MTM power spectrum also showed significant biannual phenological signals overlain on the annual cycles (e.g., Figure 2) in ~20% of pixels (largely in northern and western Eurasia; Figure 7) that were far less prevalent in the NDVI and VOD, suggesting that SIF may detect subannual photosynthetic phenologies that are not present in canopy morphology or water content of aboveground biomass.   Note that these are not mutually exclusive categories, so a given pixel can have a significant biannual phenological cycle overlain on an annual phenological cycle (e.g., Figure 2). Significance of (bi)annual phenological dynamics was defined as the observed PSD from the MTM test exceeding the 99 th percentile of PSDs from a noise-only ensemble.
The MTM power spectrum detected fewer distinct growing seasons for all three remotely sensed indices in tropical and arid regions (Figures 6 and 7). In the Amazon, VOD generally detected more distinct annual growing seasons than NDVI or SIF ( Figure 7A-C), possibly reflecting the relatively low sensitivity of VOD to the persistent cloud cover in this region [15,17,19,20]. All three indices, especially SIF and NDVI, captured clear biannual phenological dynamics in tropical forests of Africa ( Figure 7D-F). SIF and VOD also showed both annual and biannual growing seasons in parts of the Sahel, whereas NDVI did not clearly detect the biannual component of the growing season in this region ( Figure 7D-F).
In the least productive regions (phenoregions 5 and 9), SIF detected fewer distinct growing seasons than NDVI or VOD (Figures 6 and 7). This may be somewhat surprising given the presence of clearer seasonal cycles in the mean SIF phenology (Figure 3), but we suspect that this discrepancy the monthly NDVI, SIF, and VOD data. Note that these are not mutually exclusive categories, so a given pixel can have a significant biannual phenological cycle overlain on an annual phenological cycle (e.g., Figure 2). Significance of (bi)annual phenological dynamics was defined as the observed PSD from the MTM test exceeding the 99 th percentile of PSDs from a noise-only ensemble.

Recent Trends in the Spatial Extent of the Global Phenoregions
During the study period (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), the overall proportion of land area in relatively low productivity phenoregions increased at the expense of the highest productivity phenoregions ( Figure 8; Table 1). The largest losses occurred in phenoregion 4 (approximately −205,000 km 2 yr −1 ), largely consisting of temperate forests and croplands, and phenoregion 12 (−365,000 km 2 yr −1 ), largely representing tropical forests. The largest gains in global land area over the study period were in phenoregions 1 (+160,000 km 2 yr −1 ) and 2 (+89,000 km 2 yr −1 ), both of which were relatively low-productivity, high-latitude systems, and phenoregion 11 (+153,000 km 2 yr −1 ), which mostly occurred in the subtropics. These changes could represent either large-scale deforestation of the tropics [55,56] or could indicate reductions of productivity in temperate and tropical forests due to temperature-induced moisture constraints and the expansion of lower productivity ecosystems that are predominantly temperature constrained [57,58]. However, since our time series was relatively short and since we did not explicitly account for land use or land cover change, we suggest that the causes and consequences of these shifts in phenoregion area merit further research. Combining long-term NDVI and VOD datasets with a multisensor SIF synthesis [59] or with reconstructed SIF based on other optical sensors [60] could extend the length of the temporal record and allow longer-term change detection [61].
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 14 productivity, high-latitude systems, and phenoregion 11 (+153,000 km 2 yr −1 ), which mostly occurred in the subtropics. These changes could represent either large-scale deforestation of the tropics [55,56] or could indicate reductions of productivity in temperate and tropical forests due to temperatureinduced moisture constraints and the expansion of lower productivity ecosystems that are predominantly temperature constrained [57,58]. However, since our time series was relatively short and since we did not explicitly account for land use or land cover change, we suggest that the causes and consequences of these shifts in phenoregion area merit further research. Combining long-term NDVI and VOD datasets with a multisensor SIF synthesis [59] or with reconstructed SIF based on other optical sensors [60] could extend the length of the temporal record and allow longer-term change detection [61].  Table 1.

Conclusions
Determining the causes and consequences of phenological change, including its manifestations across the spectrum of ecosystem structure and function, is a key challenge in global change ecology. From the satellite perspective, defining regions with similar structural and functional phenologies (so-called "phenoregions") can improve monitoring and modeling of land surface phenology and its  Table 1.

Conclusions
Determining the causes and consequences of phenological change, including its manifestations across the spectrum of ecosystem structure and function, is a key challenge in global change ecology. From the satellite perspective, defining regions with similar structural and functional phenologies (so-called "phenoregions") can improve monitoring and modeling of land surface phenology and its change through time. Here, we demonstrated the utility of the self-organizing map method for defining a hierarchical, topologically organized set of phenoregions based on multiple satellite-based vegetation indices: NDVI (sensitive to changes in leaf area and chlorophyll content of the canopy), SIF (sensitive to photosynthetic activity), and VOD (sensitive to water content of aboveground biomass). Together, these indices capture multiple aspects of land surface phenology, allowing the delineation of phenoregions that go beyond the "greenness" paradigm [13,62] to capture a greater range of phenological variability among different ecosystems. This multisensor phenoregion approach may also provide an alternative to the land cover-based plant functional type categories that are typically used in process-based land surface models, which add considerable uncertainty to modeled estimates of carbon, water, and energy fluxes [63].
In addition to defining these phenoregions, we also highlighted the presence (and absence) of annual and biannual growing season dynamics among the phenoregions and satellite-based vegetation indices. We showed, for example, that while SIF detected fewer annual growing seasons than NDVI and VOD in many of the least productive systems, it also tended to capture more subannual growing season dynamics, suggesting that SIF may be capturing fast-varying physiological processes not observed in NDVI or VOD. We also showed that over the nine-year study period, which was limited to the period of overlap of the NDVI, SIF, and VOD data, there were significant changes in global land area of several phenoregions, notably a decline in the phenoregion representing the highest productivity ecosystems (mostly consisting of tropical forests). Combining long-term NDVI and VOD data with development of a longer and higher resolution SIF record could allow better characterization of this trend and attribution of the change to climatic versus land use drivers.