Our study shows the increasing availability of Landsat observations over the historical record affects annual summaries of vegetation greenness in seasonally snow-covered environments, leading to false greening trends that confound efforts to understand multi-decadal vegetation responses to climate change, particularly in climate-sensitive, high-alpine environments. Specifically, we confirmed that annual estimates of maximum NDVI systematically increase with the number of available observations per growing season, as previously shown for the Arctic (Berner et al. 2020, 2023), while further showing that other annual metrics (e.g., median NDVI) are similarly affected. Focusing on the European Alps, we then demonstrated how the increasing availability of observations over the past four decades led to false greening trends associated with the sampling biases and, moreover, that the magnitude of these false greening trends was highest in areas with the lowest observation density and shortest growing season length. For instance, almost 50% of the observed greening trend in high-elevation alpine grasslands was due solely to this effect, though this effect was less pronounced in lower-elevation forest ecotone communities. Annual metrics of vegetation greenness are increasingly being derived from Landsat data (e.g., Rumpf et al. 2022, Bayle et al. 2022); however, our study underscores that caution is needed when deriving and analysing annual metrics of vegetation greenness using these satellites. To rigorously quantify multi-decadal changes in vegetation greenness, and likely other biophysical parameters, it is necessary to account for the effects of changes in observation availability throughout the Landsat record.
4.1. Implications for relating greening trends to ecological processes
Arithmetically derived estimates of annual maximum NDVI are impacted by sampling bias, leading to false greening trends in subsequent time series analyses. Consequently, by not considering the effect of sampling bias, the ecological interpretation of the greening observed by Landsat turns out to be partly erroneous. For example, when compared to the actual greening observed based on our computation and those of Rumpf et al. (2022) (Figure S6), the false trends might correspond to more than 50% of the observed greening in the alpine zone (> 2400 m a.s.l.) of the European Alps. This might explain why Rumpf et al. (2022) found a significantly higher proportion of surfaces were greening (77%) compared to previous studies (55%) that used MODIS satellites (Choler et al., 2021). Nevertheless, it is crucial to note that our demonstration in no way calls into question the overall greening of the European Alps, repeatedly demonstrated by several satellites and in many different contexts (Bayle, 2020; Bayle et al., 2023; Carlson et al., 2017; Choler et al., 2021; Dentant et al., 2023; Rumpf et al., 2022), as the absolute values of the false trends we obtained are systematically lower than the trends observed without correction (Figs. 5 and 7).
Growing season length and Landsat observation availability are the two parameters that most influence the magnitude of the false greening trends (Fig. 4), both of which exhibit spatial variability (Fig. 5). In complex terrain, vegetation phenology is highly variable over short distances because of varying snow cover duration, thus the difficulty of capturing the true seasonal maximum NDVI can also vary over short distances (Choler, 2018; Dedieu et al., 2016). In addition, the number of usable observations over the growing season varies along these same gradients. During the growing season in the European Alps, cloud cover affects an average of ~ 25% of Landsat observations, but up to 50% at the highest elevations (Hu et al., 2019). Furthermore, usable observations tends to decrease with increasing elevation because of poorly performing cloud masks that erroneously assume similar temperature from all clear-sky observations in one image, resulting in misclassification of cold surfaces as clouds (Qiu et al., 2017). High-elevation alpine ecosystems are particularly affected by false greening trends because of the late snowmelt dates and limited observations, while also being among the most threatened by climate warming (Engler et al., 2011; Gottfried et al., 2012; Steinbauer et al., 2018). One of the most striking examples is that of snowbed communities, which are particularly sensitive to climate and snow cover regime changes (Hiller et al., 2005; Matteodo et al., 2016; Schöb et al., 2008) and strongly affected by false greening trends (Fig. 6). Monitoring these endangered habitats is crucial to guide conservation policies, yet long-term surveys are generally lacking, Landsat time series are an ideal candidate for tracking the long-term response of these habitats. Yet, we showed that not accounting for NDVImax sampling bias and resulting false greening trends could have detrimental consequences for understanding changes that are occurring in these habitats (Fig. 6). Beyond snowbed communities, which represent an extreme case in terms of bias magnitude, the linear increase in bias with elevation above 2400 m limit our ability to study the thermophilisation process using Landsat, as repeated plant surveys are carried out along a similar elevation gradient up to 3000 m a.s.l. (Dentant et al., 2023; Gottfried et al., 2012; Rosbakh et al., 2014).
Our study did not consider the possibility of long-term increases in growing season length or reductions in snow cover, which could increase the number of usable observations and lead to better estimates of NDVImax at the end of the Landsat time. Over recent decades, climate change has led to longer growing seasons throughout the Northern Hemisphere (Jeong et al., 2011; Park et al., 2016; Piao et al., 2019). In cold regions, many studies have highlighted the reduction of snow cover duration during recent decades with various magnitudes depending on the region and elevation (Brown et al., 2010; Callaghan et al., 2012; Hernández-Henríquez et al., 2015; Mudryk et al., 2020). For example, the Eurasian Arctic region experienced larger declines in the duration of the snow-covered period (12.6 days), compared to the North American Arctic region (6.2 days) between 1982 and 2011 (Barichivich et al., 2013), though both changes could affect the magnitude of false greening trends (Fig. 5). For global mountains, shortening of the snow cover duration is observed across elevations albeit with various magnitudes (Beaumet et al., 2021; Immerzeel et al., 2009; Klein et al., 2016; Monteiro & Morin, 2023; Notarnicola, 2022). This additional temporal and spatial variability in changing snow cover duration and growing season length are expected to add complexity that was not accounted for in our analysis.
In this study, we evaluated false and observed greening trends based on linear changes over the entire Landsat period. However, recent studies show high temporal variability in NDVI trends, with phases of greening, stability and even browning (Phoenix & Bjerke, 2016). While overall greening has been observed during the last 40 years on a global scale, there is pronounced spatial variability and local reversals toward browning seem to be increasingly observed in part due to air temperature approaching or exceeding its optimum for vegetation and increasing moisture limitation (Huang et al., 2019; Keenan & Riley, 2018; Winkler et al., 2021). In cold regions where Landsat is the most relevant tool to study these patterns (Berner & Goetz, 2022; Berner et al., 2020; Ju & Masek, 2016), local reversals from greening to browning have been documented (Phoenix & Bjerke, 2016). These observations tend to push research towards further consideration of the temporal heterogeneity and non-linear trends in vegetation greenness dynamics (Bayle et al., 2022; Vickers et al., 2016; Wolkovich et al., 2014).
The Landsat history of increasing observation is itself non-linear and occurred mostly at two moments, when Landsat 7 ETM + and Landsat 8 OLI came into operation in 1999 and 2013 respectively; Therefore, changes in the sampling bias likely occur non-linearly within the Landsat observation period. Bayle et al. (2022) studied the temporal heterogeneity of Landsat-based greening trends while not accounting for the sampling bias in Northeastern Canada and found the occurrence of two distinct waves of greening, centred around 1996 and 2011. Interestingly, both waves occurred at these Landsat transition periods, and while the first one was found to be highly correlated with rising summer temperatures, leaving little doubt to its ecological relevance, the second was not. Over their study site, \({{\Sigma }}_{{\text{O}\text{B}\text{S}}_{\text{G}\text{S}}}\) is approximately equal to 125 with GSL of 100 days, which suggests the presence of non-negligible NDVImax sampling bias (Figs. 4 and 6). Given that the trend of the second greening wave (\({{\beta }}_{\text{N}\text{D}\text{V}\text{I}\text{m}\text{a}\text{x}}\) = 0.0075 NDVI yr− 1) was close to the average bias magnitude of 0.001 NDVI yr− 1 found under these conditions in the European Alps (Fig. 4), it is possible that the result did not reflect ecological change on the ground and instead was caused solely by increasing observation availability. This specific case is an example of how if unaccounted for, this bias hampers our ability to interpret and understand ongoing ecological processes using remote sensing approaches.
4.2. False trends consideration for cold ecosystems worldwide
While our study focused on the European Alps and above-forest alpine vegetation, the combined effects of increased Landsat observation availability over time and short growing season length is found in cold ecosystems worldwide (Zhang et al. 2022). At a biogeographical scale, variation in Landsat observations comes from the combined effects of Landsat history (e.g., SCL failure) and topoclimatic patterns (see supplementary figures S7 and S8 for Arctic and Asian mountains ranges example). For instance, the North American subarctic has particularly high observation density in the 1990s compared to other adjacent subarctic regions that have similar values as observed in the European Alps (Fig. 1 and S7). Nevertheless, the North American subarctic shows contrasting patterns with higher observation frequency in the inner Canadian Plains compared to Alaska or Nunavik. This spatial variability is attributable to several local factors. For example, Alaska suffered from the data relay problem and limited tasking before 1999 (Goward et al., 2006), while two ground receiving stations were in full operation in Canada (White & Wulder, 2014), resulting in significant differences in observation availability (Ju and Masek, 2016). At the same time, continentality affects cloud cover distribution, resulting in limited observations in, for example, Nunavik (Figure S8). Similar to alpine ecosystems, growing season length varies drastically along latitude and topography gradients in sub-arctic to high arctic regions. Because the Landsat time series has increasing observation over time with large spatial variability in its magnitude, these underlying patterns could explain part of the spatial variability in Landsat-based vegetation greening trends. The AVHRR satellite time series also has inconsistent observation numbers over time (Dech et al., 2021). Hence, it is possible that differences between the Landsat and AVHRR satellites in spatial patterns of vegetation greening across North America (Ju & Masek, 2016) could partially be due to satellite-specific changes in observations availability over time. In the mountains of Asia, there are very few observations at the beginning of the time series, whereas at the end of the series the number is equivalent to that of the other cold ecosystems presented (Figure S9). Thus, the trends in the increase in the number of observations are much more acute, suggesting a potentially greater bias. Furthermore, in certain climatic contexts, such as the Austrian Alps, higher cloud cover can coincide with peak vegetation, making it more difficult to capture NDVImax than in the southern Alps. Although our analysis focused on the European Alps, the mechanisms behind the emergence of false trends (increasing number of observations over time, low observation density, and short growing season length) can be found throughout the cold ecosystems of the Northern Hemisphere (Fig. 1, S7 and S8). Overall, spatial variability in observation and phenology lead to spatial variation in the sampling bias magnitude and importance, which might hinder our understanding of spatial variability in greening trends (Myers-Smith et al., 2020).
4.3. Corrections or solutions proposed in the literature
The effects of sampling bias on phenological metrics like maximum NDVI can be addressed by either an avoidance or a correction strategy. While the sampling bias is expected to affect all terrestrial surfaces because of global patterns of increasing image availability over time (Figs. 1 and supplementary figures), we have shown that it only becomes a serious problem in certain cases (Figs. 5 and 7). For example, several studies have limited their analysis to years with sufficient observations (Carlson et al., 2017; Dentant et al., 2023; Frost et al., 2014; Raynolds et al., 2018), while others limited the study areas to overlapping tiles (Fraser et al., 2011). If skewed surfaces cannot be avoided, i.e., above 2400 m in the European Alps (Fig. 6), corrections have to be considered. Correcting for the bias means reconstructing or constructing missing information from information available elsewhere, a field in which the literature is flourishing. According to Shen et al. (2015), reconstruction approaches can be classified into three categories depending on the additional information used: (1) spatial-based methods, (2) spectral-based methods and (3) temporal-based methods. In our case, as the missing information is often outside satellite acquisitions, the first two methods are generally ineffective. Temporal reconstruction can either be based on subsequent information, assuming regular chronological fluctuation (Beck et al., 2006; Chen et al., 2004; Viovy et al., 2007), or additional information, using data fusion methods (Gao et al., 2006; Li et al., 2021; Qiu et al., 2021; X. Zhu et al., 2016). In the first case, a Landsat-based temporal interpolation approach first fits a time series model to the original Landsat NDVI data and then predicts NDVI values on regular temporal intervals (Kovalskyy & Henebry, 2012; Melaas et al., 2013; Yan & Roy, 2020; Zhang et al., 2021; Zhu et al., 2015). This approach was used by Berner et al. (2020) to correct for NDVImax estimation before computing trends in vegetation greenness across the Arctic. They developed a phenology-based approach by modelling the per pixel land surface phenology over a 17-year moving period to predict annual NDVImax and demonstrated the effectiveness of the correction. This approach has been effectively used with other spectral indices, e.g., EVI2, kNDVI; (Berner & Goetz, 2022) and is implemented in the ‘LandsatTS’ software package for R (Berner et al., 2023). In the second case, Landsat spatiotemporal fusion with another satellite of coarser spatial resolution but higher temporal resolution is a popular technique for reconstructing Landsat NDVI time-series (Chen et al., 2021; Moreno-Martinez et al., 2020). In complex landscapes only MODIS could provide adequate spatial resolution, but data is not available before 2000. AVHRR products have too coarse a resolution and have been shown to be of limited use in correcting the bias, even in a forestry context (Wang et al., 2022). With the openings of the SPOT World Heritage archives in 2015, these non-systematic acquisition products might offer interesting possibilities to densify NDVI time series before applying a temporal filter (Barrou Dumont et al., 2023; Nosavan et al., 2017).