Time-Scale Dependent Relations Between Earth Observation Based Proxies of Vegetation Productivity

Much uncertainty remains in measuring the inter-annual and longer-term dynamics of vegetation gross and net primary productivity (GPP, NPP) and the connected land carbon sink. Potential for better GPP estimation lies in newer satellite products representing different processes or vegetation states, but how they capture interannual GPP dynamics remains to be explored. Here, we differentiate shorter-and longer-term vegetation dynamics and their drivers for several Earth-observation-based vegetation proxies and a GPP estimate using time series decomposition. We find that relations between proxies can significantly differ between time scales, along land cover and climate gradients. For GPP estimated at FLUXNET sites, seasonal and multiannual slopes differ by at least 25% for half of the cases investigated, indicating considerable mismatch if multiannual relations were derived from seasonal slopes. Considering time-scale-specific sensitivities between proxies of vegetation productivity may improve estimates of interannual variability in vegetation productivity.

interannual and longer time scales. Before drawing conclusions on trends or long-term variability from seasonally dominated relations, it is important to test whether the observed relations are time-scale-invariant. To test this hypothesis, we present a systematic comparison of vegetation proxies across time scales in this study.
Different processes govern plant productivity across time scales. When using SIF as a proxy for photosynthesis, for example, studies find that it represents different physiological processes depending on the temporal and spatial scale considered (see e.g., Magney et al., 2020). At leaf to ecosystem level on daily to biweekly temporal resolution, chlorophyll fluorescence can, for example, track reduction of photosynthesis due to plant stress responses (e.g., Ač et al., 2015;Castro et al., 2020;Flexas et al., 2002;Shekhar et al., 2020;Wieneke et al., 2018, but see also Wohlfahrt et al., 2018) or can be dominated by canopy structure or absorbed photosynthetically active radiation K. Yang et al., 2018;. On seasonal scales, SIF shows a strong link to phenology-driven changes in chlorophyll content and absorbed photosynthetically active radiation (APAR, Magney et al., 2019;Wieneke et al., 2018;K. Yang et al., 2018), but can be a better proxy for photosynthetic activity than reflectance-based proxies, for example, in evergreen forest with winter dormancy Walther et al., 2016). Consequently, a relationship derived between SIF and GPP at the daily scale does not equal their seasonal relationship , and neither of them would necessarily be applicable to derive or predict the SIF-GPP relation at interannual and longer time scales. In other words, the relationship between GPP, vegetation proxies, and climate can be expected to differ substantially across time scales. A global quantitative assessment of this variation is to our knowledge missing.
A variety of vegetation proxies is available today and used to study plant phenology, physiology, and their relation to GPP. When comparing such proxies, it is important to note that they represent different aspects of plant physiology. For example, GPP FC and SIF approximate photosynthesis-related fluxes of carbon and energy (Badgley et al., 2017;Frankenberg & Berry, 2018), while normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), or leaf area index (LAI) represent vegetation states related to greenness or leaf area of standing vegetation, and hence photosynthetic capacity (Huete, 1997;Tucker & Sellers, 1986). Further differences exist by design, where, for example, NDVI and EVI are bounded variables, while SIF is not confined by an upper boundary. NIR V can be calculated based on radiance, in which case it behaves more like a flux variable, or reflectance, in which case it is bounded and unitless like EVI or NDVI. These differences lead to known nonlinear relationships between proxies: One prominent example is the nonlinear GPP-NDVI relation over evergreen forest with winter dormancy (e.g., Magney et al., 2019), as well as at high GPP values, caused by saturation of NDVI over densely vegetated areas (Frankenberg et al., 2011).
In more general terms, satellite-based proxies of vegetation capture different aspects of plant physiology and phenology depending on their design, as well as their temporal and spatial scale of observation. Thus, differentiating the relation of satellite-based vegetation proxies between shorter and longer temporal scales could add to better understanding the evolution of GPP dynamics in the future. In a previous study, we characterized major temporal modes of vegetation variability in NDVI time series, and found that correlations between NDVI and climate variables could even invert between time scales (Linscheid et al., 2020). We here extend this time scale analysis to several newer vegetation proxies to understand their relations across time scales, and how these change with climate and land cover. We decompose time series by Discrete Fast Fourier Transform (FFT) and extract subsignals at subseasonal, seasonal and multiannual time scales. We investigate when and how linear relations between vegetation proxies differ between subseasonal, seasonal and multiannual time scales globally: We determine time-scale specific slopes between GOME-2 SIF, FLUXCOM GPP (GPP FC ), MODIS NIR V , GIMMS NDVI, and MODIS EVI. We then test how slopes vary in relation to time scale, land cover and climate to infer when and where scale-specific slopes should be considered for better understanding dynamics of vegetation productivity. We corroborate our findings for site-scale GPP estimates from 14 FLUXNET sites and with an ensemble of SIF products.

Global Datastes
We compared global datasets of satellite-based vegetation proxies. Of the data-driven FLUXCOM RS GPP product (GPP FC ) we used the full ensemble mean monthly data (Jung et al., 2020;Tramontana et al., 2016). Throughout the text, we explicitly differentiate between this data set approximating global GPP (abbreviated as GPP FC ) and instances where we refer in more general terms to actual/true gross primary productivity (abbreviated as GPP). For sun-induced fluorescence (SIF), we used data from GOME-2 SIF retrievals originally at monthly temporal and 0.5° spatial resolution (GFZ, Köhler et al., 2015). Enhanced vegetation index (EVI, Huete, 1997) from the BRDF-corrected MODIS MCD43C4.006 reflectances were used (Didan, 2015), originally at daily temporal and 0.05° spatial resolution. Near-infrared reflectance of vegetation (NIR V ) was also calculated based on BRDF-corrected MODIS MCD43C4.006 reflectances after Badgley et al. (2017). Normalized Difference Vegetation Index (NDVI) was retrieved from GIMMS 3g v1 (NCAR, 2018;Pinzon & Tucker, 2014;Tucker & Sellers, 1986;Tucker et al., 2005), originally at biweekly temporal and 1/12° spatial resolution. The overlapping period of all products were the years 2007-2015. For comparability with climate datasets, the vegetation data were resampled by weighted averaging to monthly temporal and 0.5° spatial resolution where necessary.
We additionally employed an ensemble of SIF datasets including NASA GOME-2 SIF MetOp-A v28 level3 (Joiner et al., 2014), the KNMI SIFTER L2 v2 product van Schaik et al., 2020) and a calculation of a SIF daily integral from instantaneous GFZ GOME-2 SIF mentioned above. These estimates of daily SIF were derived from the instantaneous soundings using the daily integral of the cosine of the sun zenith angle scaled by the cosine of the sun angle at the observation time (Frankenberg et al., 2011).
We investigated global patterns of time-scale specific slopes in relation to land cover and climate, as well as along aridity and temperature gradients. To this end, we used several categorical maps to split the data set. If not already at the correct resolution, datasets were aggregated to 0.5° spatial and, if applicable, monthly temporal resolution. A map of Koeppen-Geiger climate zones (Kottek et al., 2006) was obtained through splitting by first letter of the classification (A, equatorial, B, arid, C, warm temperate, D, snow), or additionally by second letter to account for precipitation seasonality (W, S, f, s, w, m). A map of land cover type was calculated from the synergistic land cover product SYNMAP land cover fractions (Jung et al., 2006). To retain only the most homogenous pixels for each land cover type, data below 70th quantile and below 40% land cover fraction was discarded for each major land cover type. We also used SynMap to calculate the fraction of vegetation cover and fraction of tree cover per pixel. As temperature data set, we used ERAinterim v2 air temperature (Dee et al., 2011) and calculated mean temperature over the whole time period. The aridity index, a measure of relative water-or energy limitation of a system, was calculated as the ratio of surface net radiation (SNR) and total precipitation (TP). It was calculated from ERA5 TP and mean daily SNR for each time step via the formula aridity = 30*SNR/2.45*TP (Allen et al., 1998). Subsequently data were temporally aggregated by median over the whole time period. For all plots showing results summarized by category or binning, data were only displayed if at least 30 pixels remained in the respective bin.

Datasets for FLUXNET Site-Scale Analysis
For an analysis at FLUXNET sites, site-scale estimates of GPP from eddy covariance came from the FLUXNET 2015 collection, with standard ONEFlux processing (Pastorello et al., 2020). GPP estimates were derived via the daytime NEE partitioning method  with standard MDS gap filling and quality control (Pastorello et al., 2020). Sites were selected to have long time series, which maximized eddy covariance and MODIS data availability as described below.
We extracted BRDF-corrected MODIS surface reflectances from the MCD43A4 product at 500m pixel size and daily sampling (Schaaf & Wang, 2015b). The accompanying MCD43A2 product (Schaaf & Wang, 2015a) was used for quality control. Both data sets were reprojected to a regular grid at the tower locations. Quality checks reject data points with the BRDF_Albedo_Band_Quality_Band1-7 > 2 and BRDF_Albedo_LandWatertype > 2 and allow only observations where the flag Snow_BRDF_Albedo indicates snowfree conditions. From these quality filtered surface reflectances, spectral indices EVI, NDVI, and NIRv are computed for the pixel containing a tower.
MODIS and GPP datasets were aggregated to monthly resolution and processed analogously to global datasets (see below for details). We identified 14 Fluxnet sites running over the whole MODIS period, where enough MODIS data was available in the pixel containing the tower to represent the inherent seasonal cycle and interannual variability and ensure proper representation of seasonal and multiannual dynamics, and where gapfilling did not introduce artificial interannual variability. These sites were: AU-Tum, BE-Vie, CA-TP4, DK-Sor, FI-Hyy, FR-Pue, IT-Col, RU-Fyo, US-MMS, US-Me2, US-PFa, US-Ton, US-UMB, and US-Var.

Decomposition and Regression
All analyses were performed within the Earth System Data Lab framework (ESDL.jl, Mahecha et al., 2020) using the Julia programming language (Bezanson et al., 2017). Time series were gapfilled per pixel in two steps: first using a fourth order polynomial, subsequently filling remaining gaps at the ends of each time series with values from the mean annual cycle. Most gaps occurred in northern latitudes ( Figure S1). Time series were decomposed based on Discrete Fast Fourier transform (FFT, Brockwell & Davis, 2006), followed by binning of frequencies into subseasonal, seasonal, and multiannual bands as previously described (Linscheid et al., 2020). Briefly, for each pixel the mean and linear trend are removed from the time series. The remaining signal is decomposed by FFT and binned as follows: seasonal signals are reconstructed from frequencies of 0.9-1.1 years and two seasonal harmonics (0.5 and 0.33 years), the subseasonal (or short term) bands contain the remaining frequencies <0.9 years, and the multiannual band contains all frequencies >1.1 years (interannual and longer scales). When using FFT, the first and second harmonic of the seasonal cycle (0.5 and 0.33 years) carry residual information necessary to recreate a good representation of a seasonal cycle and are hence included in the seasonal band here. All remaining (i.e., not overlapping) frequencies <0.9 years are then attributed to the subseasonal signal. Each subsignal still contains the same number of data points as the original time series, that is, time scales are not attained by data aggregation, but separated by partitioning of different frequency bands. Summing all extracted subsignals, trend, and mean for each point in time returns the original signal. FFT could only be performed on pixels without missing data points. Gapfilled data points were masked following the FFT decomposition step to base subsequent analyses only on original data points.
To quantify time-scale specific relations between variables, we performed pixelwise standard major axis regression (SMA, also known as total least squares regression, TLS), for each pairwise combination of variables at each of the extracted time scales. SMA is needed given that both variables are expected to contain errors and the method minimizes errors in both x and y direction (Legendre & Legendre, 2012). Thus, the slope is hardly affected by the correlation: The SMA slopes are identical then to the first principal component of a bivariate data set. In SMA, slopes do not depend on which variable is the predictor, that is, slope SMA (GPP vs. SIF) = 1/slope SMA (SIF vs. GPP). For each pairwise comparison of vegetation proxies, the naming convention follows the order y-axis versus x-axis, that is, the first variable is plotted on the y-axis, the second on the x-axis (e.g., SIF-GPP FC means GPP FC is plotted on the x-axis, SIF on the y-axis). With GPP FC plotted on the x-axis, high slopes indicate high sensitivity of the respective proxy to changes in GPP FC (a change in y with change in x), while slopes close to zero indicate low sensitivity (little change of y with change in x).
For each variable pair, the numerical value of their slope depends on the variable units and order of magnitude, and thus cannot be directly compared across pairs. To compare slopes of different variable pairs, and identify regions with varying slopes across time scales, we thus computed the ratio of the slopes of multi-annual signals to the slope of seasonal signals at each pixel (slope ma /slope seas ). These slope ratios are unitless and thus comparable across variable pairs, and indicate where slopes are higher at seasonal or multiannual time scales. The slope ratios illustrate the relative partitioning of variance between seasonal and multiannual signals (for illustration see Figure S2). For convenient interpretation, we summarize the analysis as log 10 (abs (slope ma /slope seas )), so that logarithmic slope ratios below zero (blue in the figures) indicate that slopes are higher in the seasonal than multiannual band: This originates from relatively less multiannual variability captured in the y-variable than in the x-variable; that is, the "sensitivity" of y to x is higher on seasonal scale ( Figure S2a). Conversely, positive slope ratios (red in the figures) indicate higher multi-annual slopes and thus relatively higher multiannual than seasonal sensitivity of y to x ( Figure S2b).

Data Filtering
We applied strict filtering to mask pixels that could cause spurious results. Signals with very small amplitude (low spectral power in this time scale) could cause aberrant slopes or slope ratios. For each variable, we thus masked all pixels where the amplitude of the decomposed signal accounted for less than 5% of the original signal. For the same reason, we further removed pixels accounting for the 4% most extreme amplitude ratios for each variable pair (retaining quantiles 0.02-0.98 of their global distribution). SMA regression should only be computed when a significant correlation between variables exists (Legendre & Legendre, 2012); we therefore masked all pixels where the correlation significance showed a p-value > 0.1, where correlation significance was computed by One-Sample Z-Test. For comparison across variable pairs, we then combined the pairwise masks into a common mask, only retaining pixels that passed for all variable pairs and all time scales. This mask retained 26% of vegetated land pixels ( Figure S3).

Relations Between EO Proxies Change Across Time Scales
Time series of NDVI, EVI, NIR V , SIF, and GPP FC were decomposed at each pixel to extract modes representing subseasonal, seasonal, and multiannual time scales of variability. As illustrated in Figure 1, these subsignals all contain a certain variability on their respective scale. Spectral approaches have been shown to extract major modes of variability from ecological time series (Katul et al., 2001;Mahecha et al., 2007;Stoy et al., 2009) and allow analysis of scale-specific variability by, for example, removing the effects of confounding factors varying on a seasonal time scale from an instantaneous relationship (Mahecha, Reichstein, Carvalhais, et al., 2010).
We first illustrate the calculation of scale-specific slopes in a comparison between SIF and GPP FC (Figure 1, Figure S4), where we find global patterns highlighting differential slopes across time scales. With GPP FC on the x-axis as chosen here, higher slope values indicate a larger change (higher "sensitivity") of the y-variable (here: SIF) to changes in GPP FC . Overall, spatial patterns of slope intensity differ between the seasonal and multiannual time scale (Figure 1c). Except for south-western Australia, slopes are positive, and highest in northern North America, northern Eurasia, Europe, as well as subtropical Africa and South America.
Much uncertainty in prediction of GPP dynamics is associated with interannual variability of GPP. Some of this uncertainty may be alleviated by better characterizing and constraining predictions of multiannual GPP dynamics with observations. We thus here choose to focus on vegetation variability on multiannual time scales, and how and where multiannual relations between proxies differ from seasonal relations. In order to compare slopes The logarithmic slope ratio of multiannual over seasonal slopes (bottom). Ratios of slopes are shown as log10 (absolute [slope ma /slope seas ]), that is, red means higher (absolute) slope in the multiannual band, blue means higher (absolute) slope in the seasonal band. Logarithmic slope ratios below zero (blue in the figures) indicate that slopes are higher in the seasonal than multiannual band. This originates from relatively less multiannual variability captured in the y-variable than in the x-variable; that is, the sensitivity of SIF to GPP FC is higher on seasonal scale. Conversely, positive slope ratios (red in the figures) indicate higher multiannual slopes and thus relatively higher multi-annual than seasonal sensitivity of SIF to GPP FC (see also Figure S2). Slope units are (mW sr −1 nm −1 /gC d −1 ); gC: gram Carbon.
between time scales, we calculated the ratio between multiannual and seasonal slopes for each grid cell (Figure 1d), indicating where slopes are higher at the seasonal or multiannual time scale. In drier regions of the globe, seasonal slopes between SIF and GPP FC are generally higher (or similar to) multiannual slopes (blue regions), while multiannual slopes are higher in wetter regions (red regions). Overall, we find that slopes between SIF and GPP FC vary with time scale, and their ratios show spatial patterns related to climate and land cover as elaborated in the next paragraph ( Figure 2).
Specifically, slopes between SIF and GPP FC vary with land cover type, as well as along gradients of temperature and aridity. Slope patterns along tree cover and vegetation cover are similar across scales (though different in magnitude, Figure 2a), while marked differences exist for climate and land cover splits across time scales (Figures 2b and 2c). The largest slopes are found between subseasonal signals, especially for the more arid, as well as warm and relatively humid climates, and nontree vegetation (Figures 2a-2c), indicating a higher variability of SIF, or lower variability of GPP FC , in these regions and thus a higher sensitivity of SIF to GPP FC . A marked exception are evergreen trees in KG class A (Figure 2c), but note that most pixels in this class belong to seasonally (d-f) same splits as in (a-c). but displaying logarithmic slope ratios of multiannual over seasonal scale. Colors represent area-weighted mean values per class of pixels. Slope ratios are shown as log10 (absolute [slope ma /slope seas ]), that is, red means higher (absolute) slope in the multiannual band, blue means higher (absolute) slope in the seasonal band. At least 30 nonmissing values were required per bin, otherwise bins were masked (indicated by hatching). A, equatorial climate; B, arid climate; C, warm temperate climate; D, snow climate; Tree_dec, deciduous tree cover; Tree_evg, evergreen tree cover. dry ecosystems (As, Aw), while fully humid tropics could not be assessed due to too few pixels remaining after filtering. Chen et al. (2020) previously identified decreasing GPP/SIF ratios from cold-and-wet to hot-and-dry climate, concluding that moisture availability critically influences the GPP/SIF ratio via effects on stomatal conductance. Our results corroborate this pattern for the subseasonal scale, but not clearly for longer time scales: Seasonal and multiannual slopes are generally lower and more uniform than short-term slopes. Seasonal and multiannual slopes are additionally high in temperate, rather arid climate. SIF thus appears generally as a more sensitive proxy for GPP FC in more arid, nontree ecosystems, as well as tropical evergreen forests that experience a dry season. SIF may have a higher inherent noise level than other vegetation proxies, thus short-term relations could be diverging due to different noise levels especially in low-productivity ecosystems. We are however using SMA to compute slopes, which is less sensitive to noise than regular linear regression. Additionally, we see high sensitivities between SIF and GPP FC also for evergreen trees and crops, two potentially very productive ecosystem types. In complex canopies, a decrease in canopy escape ratio of fluorescence could additionally cause lower observed SIF relative to observed GPP Guan et al., 2016), and thus lower slopes.
Ecosystem-specific slopes between SIF and GPP FC have been reported previously between different land cover types Guanter et al., 2012;Parazoo et al., 2014;Wood et al., 2017) or covarying with climate . On the other hand, within the same land cover type, previous studies at eddy covariance sites have found consistent slopes between GPP FC and SIF for hourly, daily, weekly and up to quarterly aggregation windows Magney et al., 2019;Turner et al., 2020;Wood et al., 2017). At the longer temporal scales assessed here, we do not generally find this time scale consistency: Seasonal slopes are higher than multiannual slopes for nontree ecosystems of all climate types (Figures 2d and 2f, blue colors) and generally in more arid climates (Figures 2e and 2f), while either higher or similar for crops, grasses and shrubs (Figure 2f, red and white areas). In all other areas, multiannual slopes are higher than seasonal slopes, indicating a higher sensitivity of SIF to GPP FC on multiannual than seasonal scale. This may be interpreted as SIF capturing relatively more multiannual variability than GPP FC in semiarid and nontree ecosystems. Note, however, potential shortcomings in long-term SIF retrievals and GPP FC IAV: Long-term SIF dynamics should be interpreted with caution, as trends derived from SIF signals have been attributed to sensor drift rather than vegetation dynamics (Zhang et al., 2018, but note we do not assess trends here). Several other known issues exist, including a swath width adjustment for GOME-2 in July 2013 which may affect analyses of the full-length time series (Parazoo et al., 2019;van Schaik et al., 2020). To minimize effects of biases introduced by the SIF retrieval, we repeated the same analysis with an ensemble of SIF products and found similar results across all of them ( Figure S5).
Similarly, interannual variability of GPP in FLUXCOM RS is likely underestimated, and more uncertain than seasonal fluxes (Jung et al., 2020;Tramontana et al., 2016). An underestimated multiannual GPP FC signal would impact the slopes and slope ratios calculated here: With smaller than expected changes in GPP FC (on the x-axis), for equal change in the y variable, multiannual slopes would be artificially increased. Additionally, differences across land cover classes may be scale-dependent and/or in part be due to biases in GPP products (Sun et al., 2018;Turner et al., 2020). We thus expand our analysis to several other vegetation proxies in the next section and compare the above results to an analysis of site-level FluxNet GPP estimates.

Comparison of Time-Scale Specific Relations Across Vegetation Proxies
To test whether the conclusions drawn above are generalizable using other vegetation proxies, we further compared EVI, NIR V and NDVI as observed satellite proxies against each other, and against SIF and GPP FC (Figures S6-S14). For SIF versus GPP FC , we found that slopes were highest on subseasonal scales. In comparison, relations of NDVI, NIR V , and EVI with GPP FC show the highest slopes in the multiannual band. NIR V , EVI and NDVI are slightly more sensitive to GPP FC on subseasonal than seasonal scales which is somewhat unexpected, because GPP would generally vary more on shorter time scales (hourly-weekly scales, e.g., caused by fluctuating land-surface temperature and radiation), while EVI and NDVI vary more slowly. Since the datasets are at monthly resolution, this expectation raised by studies operating on subdaily to weekly scale may simply not be visible here. Additionally, SIF and GPP FC carry information on actual photosynthetic activity, while the vegetation indices mostly carry information on phenology. A multiplication of EVI, NIR V , and NDVI with photosynthetically active radiation (PAR) would produce estimates of photosynthetic activity more similar to SIF or GPP, and could change their variance captured especially on the seasonal scale, which could also influence their scale-specific slopes toward SIF and GPP FC . NIR V and EVI show the most homogeneous slopes with GPP FC along gradients of land cover and tree cover fraction ( Figure S15), without the previously observed, systematically higher slopes for pixels with tree cover <10%. This may encourage generalizability of a NIR V -GPP FC relationship across time scales, yet this homogeneity was diminished when partitioning by climate and land cover class (Figures S16-S17). Overall, our results underscore that proxies such as SIF representing ecosystem fluxes capture short-term (subseasonal) variability well, while state variables are more representative of productivity variability on longer time scales. Additionally, SIF shows similar agreement with GPP FC on the longer time scales as EVI and NIR V .
For slope ratios across the variable pairs (Figure 3, Figures S18-S28), a gradient emerges: Slope ratios are highest with GPP FC on the x-axis (Figure 3, left column) due to GPP FC generally capturing less multiannual relative to seasonal variability than the other vegetation proxies, causing higher multiannual than seasonal slopes. This Figure 3. Pairwise slope ratios for all combinations of variables split by fraction of vegetation cover and fraction of tree cover. The naming convention follows "y-variable versus x-variable." Ratios of slopes are shown as log10 (absolute [slope ma /slope seas ]), that is, red means higher (absolute) slope in the multiannual band, blue means higher (absolute) slope in the seasonal band. At least 30 nonmissing values were required per bin, otherwise bins were masked (indicated by hatching). EVI, enhanced vegetation index; GPP FC , FLUXCOM gross primary productivity estimate; NDVI, normalized difference vegetation index; NIR v , near-infrared reflectance of vegetation; SIF, sun-induced fluorescence. pattern is stronger with NDVI than previously seen for SIF, but similar for EVI and NIR V . The difference between tree and nontree vegetation is most pronounced in the NDVI-GPP FC relationship, possibly related to saturation effects of NDVI at high GPP in forest ecosystems, which could dampen seasonal slopes and thus create higher slope ratios at higher tree cover fraction. Similarly, a less pronounced seasonal cycle of NDVI than GPP FC in evergreens would also have this effect.
From left to right in Figure 3, slope ratios decrease with GPP FC , SIF, EVI, and NDVI as x-variable, linked to a decreasing amount of multiannual relative to seasonal variability captured with NDVI > EVI, NIR V , SIF > GP-P FC (Figure 3, Figure S28). Despite the general shift toward higher multiannual slopes with GPP FC as x-variable, global patterns along gradients of aridity, temperature and vegetation type were similar across pairs of vegetation proxies, with exception of NIR V versus EVI, two datasets previously shown to be closely linked (e.g., Doughty et al., 2021). Our analysis thus indicates that differences in subseasonal, seasonal, and multiannual relations exist across satellite-based proxies of plant productivity. The overall finding that relations between vegetation proxies differ across time scales emphasizes the potential of differentiating multiannual and seasonal vegetation dynamics in understanding long-term vegetation growth and productivity.

Comparison With FLUXNET Site-Scale GPP
The analyses outlined above pointed to scale-specific slopes between vegetation proxies (Figure 3), but cannot rule out an effect of biased GPP FC interannual variability on the scale-specific slope comparison (Figures 1  and 2). Additionally, FLUXCOM GPP is an upscaled product and comparison to site-level GPP estimates may be advantageous. We thus additionally conducted an analysis at 14 long-running FLUXNET sites where more direct measurements of ecosystem carbon fluxes are available from the start of the MODIS data period (see Figure S29 for an example). In this data set, multiannual differed from seasonal slopes in about half of the cases ( Figure S30). We compared these results to our global results by extracting data for each nearest 0.5° pixel from the global map ( Figures S31-S32). Slope ratios differed between site-and pixel-level, with pixel-level slope ratios showing higher scale-specificity of slopes than site-level data ( Figure S31). The site-level results corroborate that differences can exist in the relation between GPP and different vegetation proxies across time scales, albeit suggesting a lower magnitude of the effect than the results above based on global gridded datasets. They also highlight, that in particular the aggregated FLUXCOM product misses out on variability on longer time-scale compared to site level observation ( Figure S31). To what extent this is intrinsic to the machine learning ensemble approach or caused by the spatial aggregation to 0.5° will need further study.

Limitations and Outlook
This study investigates broad-scale global relations between satellite-based proxies of vegetation productivity, with an emphasis on seasonal and longer time scales. The coarse spatial and temporal resolution of the datasets does not allow us to capture processes below monthly time scales or below ecosystem spatial scale. Based on the overlapping data period of eight years, the multiannual relations captured here represent multiannual, but not multidecadal dynamics in plant productivity. Due to the short time series, a finer spectral separation especially of long-term modes of variability was not meaningful, but could with longer time series provide additional insight into covariation of vegetation with climatic modes of multiannual to decadal variability. Furthermore, vegetation proxies derived from MODIS are not independent of each other, which may affect our results. However, a number of different, independent data sources used here overall point to similar conclusions: for example, analyses based on FLUXNET GPP partially supported our results at local scale based on a MODIS-independent potential GPP estimate.

Global Coverage
Most of the tropics are masked in our analysis due to constraints set by the data itself. Data coverage in tropical areas is often sparse, a general issue seen for optical remote sensing due to factors such as cloud cover impacting satellite retrievals. In our case, this affects especially the all-vs-all comparison including EVI, NIR V , and NDVI ( Figure 3). One advantage of SIF retrievals is their reduced sensitivity to cloud cover, which allows better coverage of measurements over tropical forests. In our analysis of SIF versus GPP FC (Figures 1c and 1d), more pixels in tropical climate are retained and do also show scale-specific slope ratios (Figure 1d). Hence, the effect observed elsewhere is probably also relevant in the tropics, although its quantification would require a more detailed analysis and employ a different set of data. Efforts to provide harmonized, analysis-ready high-resolution data cubes for tropical regions are under way and form an important basis for such analyses (Estupinan-Suarez et al., 2021), for example, strengthening our knowledge on seasonality of tropical vegetation (Estupinan-Suarez et al., 2021). With the era of extended radar-based biomass measurements and continuous extension of TROPOMI SIF time series, we hope that more analyses on high quality data retrieved over the tropics will become possible over the next years.

Data and Method Limitations
We employ linear statistical methods to allow comparison across time scales based on slope ratios. SMA regression is used since it is less affected by noise than regular linear regression and minimizes errors in both x and y direction. In order to minimize spurious results, we further apply strict filtering criteria like a correlation significance filter, amplitude filter, and discarding the most extreme slope ratios. Despite comparing several ensembles of satellite-based vegetation proxies and different sources of GPP estimates, uncertainties remain: Slopes are more uncertain in areas of low ecosystem productivity and at multi-annual scales as exemplified by their higher confidence intervals ( Figures S4, S6-S14). We conducted significance tests to determine whether medians of slope ratio distributions were significantly different from zero and found a consistent general shift between proxies (Figures S18-S27). Despite stringent filtering applied, it can be difficult to fully tell physiological effects apart from data characteristics. Overall, observational capacities are not yet at a stage of allowing us to draw uniformly clear quantitative conclusions on scale-specific relations of vegetation proxies across temporal scales. Still the different analyses all suggest the importance of scale-specific sensitivities in understanding vegetation productivity over time.
Processes on different time scales can interact in, for example, additive or multiplicative ways (Wu et al., 2008). A prominent example is data heteroskedasticity across seasons, where the variance among deseasonalized daily GPP values is higher during the main growing period compared to other times of the year (Lasslop et al., 2008;Menzer et al., 2013); in that case, the seasonal time scale affects the variance of the daily time scale. With regard to different vegetation proxies, one could assume that data for GPP and other fluxes contain multiplicative effects such as the one described, while more slowly changing state variables such as NDVI would not show this behavior. This still needs to be assessed for the relation between seasonal and multiannual bands. We did not assess the interaction between time scales here, but such analyses may hold additional useful information for time-scale-resolved understanding of GPP and should be pursued in the future. More detailed analyses on uncertainties in the decomposition technique would be necessary, for example, in the form of Iterative Amplitude Adjusted Fourier Transform (IAAFFT, Mahecha, Reichstein, Jung, et al., 2010;Schreiber & Schmitz, 1996) to test linearity or nonlinearity of the time series and the uncertainty of its retrieval.

Implications for Long-Term Vegetation Dynamics
State-of-the-art terrestrial biosphere models (TBMs) are often calibrated and evaluated on complete time series, which tend to be dominated by dynamics of the seasonal cycle. This may compromise their ability to capture short term, as well as interannual, and longer relations (Mahecha, Reichstein, Jung, et al., 2010). While most variables investigated here are not directly calculated within TBMs, our findings should still be valid as vegetation indices such as NDVI strongly correlate with fPAR or LAI, and can be calculated from fPAR model outputs via radiative transfer models (RTMs). Direct calculation of SIF within land surface models has also been performed (e.g., Thum et al., 2017), but is still not at a standard operational level in most TBMs. To improve representation of inter-annual vegetation variability in TBMs, they should be analyzed for their representation of longer-term GPP dynamics, and whether potentially seasonally dominated model parameterization are valid for estimating long-term carbon balances (Braswell et al., 2005;Mahecha, Reichstein, Jung, et al., 2010). Recent results of spectral Granger causality analysis based on CMIP5 models for example, indicate a partial over-estimation of climate controls on vegetation dynamics at inter-annual scales in the models compared to observations (Claessen et al., 2019).
Newer methodology also holds promise for improved modeling of dynamic systems with multiscale properties. Seen from a multiple linear regression perspective, differentiated relations at inter-and multiannual time scales would need to explicitly be addressed in model training, optimization and evaluation steps (Mahecha, Reichstein, Jung, et al., 2010) to test how much predictability of interannual variability in vegetation productivity can be improved this way versus addressing missing long-term drivers. Machine-learning-based methods can implicitly learn relationships at different scales or frequencies (Cao et al., 2020;Mendez-Jimenez & Cardenas-Montes, 2018) and may be a promising complementary approach to address time-scale dependencies in vegetation productivity dynamics and their drivers.

Conclusion
Previous studies have found constant linear scaling between vegetation proxies at time scales from hours to several months. We investigate this relationship at seasonal to multiannual scale and find that relations between vegetation proxies can diverge between these longer time scales. We find time-scale specific slopes between variables with different magnitude, but consistently for the ensemble of vegetation proxies tested here. Across satellite-based proxies, ratios of multi-annual to seasonal slopes vary between vegetation types, and along gradients of tree cover and climate. We find a gradient of multiannual relative to seasonal variability captured with NDVI > EVI, NIR V , SIF > GPP FC at the spatial and temporal resolution tested here. State-of-the-art approaches often investigate integral time series where diurnal or seasonal cycles dominate. Our study indicates that seasonal relationships do not generally hold at shorter and interannual scales. Thus, understanding and prediction of longer-term dynamics in GPP may need to explicitly take into account the multiannual signal components and their relations to Earth observation based predictors.