Greening vs browning? Surface water cover mediates how tundra and boreal ecosystems respond to climate warming

Climate warming in northern high latitudes has progressed twice as fast as the global average, leading to prominent but puzzling changes in vegetation structure and functioning of tundra and boreal ecosystems. While some regions are becoming greener, others have lost or shifted vegetation condition as indicated by a browning signal. The mechanisms underlying this ‘greening or browning enigma’ remain poorly understood. Here we use multi-sourced time-series of satellite-derived vegetation indices to reveal that spectral greening is associated with reductions in surface water cover (i.e. fraction of surface water bodies), whereas spectral browning is linked to increases in surface water cover. These patterns are consistently observed from both 30 m resolution Landsat data and 250 m resolution MODIS data on the basis of grid cells sized of 1, 2 and 4 km. Our study provides, to our knowledge, the first biome-scale demonstration that interactions between vegetation condition and water cover change can explain the contrasting trajectories of ecosystem dynamics across the northern high latitudes in response to climate warming. These divergent trajectories we identified have major implications for ecosystem functioning, carbon sequestration and feedbacks to the climate system. Further unraveling the interaction between vegetation and surface water will be essential if we are to understand the fate of tundra and boreal biomes in a warming climate.


Introduction
Northern high latitudes are warming twice as fast as the global average (Pachauri et al 2014, Post et al 2019. This warming has been linked to changes in permafrost condition (Smith et al 2005), hydrology (Walter et al 2006), vegetation (Myneni et al 1997, Piao et al 2020, and disturbance regimes (Guindon et al 2018). Disentangling the complex interactions between these abiotic and biotic changes is a fundamental step to anticipate the long-term dynamics of tundra and boreal ecosystems.
Satellite-derived spectral greening (i.e. positive trend in vegetation indices (VIs)) has been one of the most conspicuous changes associated with climate warming in northern high latitudes (Myneni et al 1997, Berner et al 2020, Myers-Smith et al 2020, Piao et al 2020. This long-term greening trend has been interpreted as enhanced vegetation growth driven by temperature rise and correlated with changes in vegetation cover (Elmendorf et al 2012), biomass (Forbes et al 2010), productivity (Berner et al 2020), abundance and height (Elmendorf et al 2012), as well as prolonged plant growing season (Zeng et al 2013, Park et al 2016. Comparable trends have been documented by field observations. For example, increased abundance and height of vascular plants have been recorded across permanent plots in the Arctic tundra (Elmendorf et al 2012). Also long-term field monitoring has demonstrated increases in shrub biomass, cover, height and abundance, as well as advanced spring phenology, and a northern expansion of shrub distribution ranges in Siberian, Alaskan and northern Canadian tundra (Frost and Epstein 2014). In the boreal biome, spectral greening has been attributed to increased forest productivity (Beck et al 2011), to transitions from coniferous to deciduous forests, as well as to a northward shift in the distribution range of trees (Wang and Friedl 2019).
Despite the widespread greening trends in tundra and boreal ecosystems, there is also evidence showing opposite trends of vegetation change inferred from satellite-derived spectral browning (i.e. negative trend in VIs). Although in some cases spectral browning could be at least partly caused by sensor degradation (Guay et al 2014), there is much field evidence of vegetation declines across certain areas of tundra and boreal ecosystems (Phoenix and Bjerke 2016). Tree mortality caused by more frequent and intensified wildfires and insect outbreaks explains some of these changes (Gamon et al 2013, Bjerke et al 2017, National Academies Of Sciences 2019, Wang and Friedl 2019. However, in areas where these disturbances are apparently absent, spectral browning and vegetation decline have been related to winter warming, frost drought and waterlogged conditions by thermokarst development Bjerke 2016, National Academies Of Sciences 2019).
Although vegetation may recover from extreme winter warming (Bokhorst et al 2011) and fire (Bret-Harte et al 2013) perturbations, it has been hypothesized that potential positive feedback mechanisms may lead to abrupt and persistent changes in vegetation structure in boreal and tundra ecosystems (Scheffer et al 2012). In tundra ecosystems, for example, a well-known feedback involves snow-shrub interactions. The presence of shrubs can augment snow accumulation in winter, raising winter soil temperature and facilitating shrub survival (Sturm et al 2001). A substantial loss of shrub cover may therefore disrupt this self-maintenance mechanism and new feedbacks may develop. In both Siberian and Alaskan sites, thawing of ice-rich permafrost resulting from rising summer temperatures can lead to waterlogged conditions resulting in shrub and tree mortality (Osterkamp et al 2000, Karlsson et al 2011, Frost and Epstein 2014. Once aboveground woody vegetation decreases, positive feedbacks could potentially unfold to further increase permafrost thawing. Small-scale in situ experiments have demonstrated that removal of shrub cover can result in summer thawing of the top permafrost and, in turn, leading to soil subsidence and increases of surface water cover which further facilitate permafrost thawing and prevent the reestablishment of shrubs (Blok et al 2010, Nauta et al 2015. Although the abovementioned lines of evidence suggest that surface water cover may play an important role in the complex interactions that determine the direction of vegetation change, there has been yet no global assessment of how changes in surface water cover, vegetation condition and climate interact with each other across northern high latitudes. In this study, we address this gap by examining the relationships between trends in surface water cover and spectral greening and browning across the tundra and boreal biomes during 2000-2015. Building on the results from small-scale in situ experiments in tundra (Blok et al 2010, Nauta et al 2015, we expected that negative associations between the trends in surface water cover and spectral greenness could be detectable at larger spatial scales, and therefore hypothesized that greening trends across high latitude landscapes were likely associated to reductions in surface water cover whereas browning vegetation trends were associated to wetter conditions reflected by increases in surface water cover. For surface water cover, we focused on the inland surface water bodies (e.g. ponds), that can be consistently detected through Landsat images at a resolution of 30 m (Pekel et al 2016).

Methods
We explore how changes in surface water cover relate with trends in spectral greening and browning (i.e. linear changes in the annual mean of VIs during the growing season) across the tundra and boreal biomes. Using time-series of satellite-derived annual surface water cover (Pekel et al 2016) and several VIs (including the normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI), derived from both MODIS and Landsat data), we calculated the temporal trends of spectral greenness and surface water cover between 2000 and 2015 on the basis of grid cells of varying sizes. We used linear regression modeling to examine the bivariate associations between the resulting grid-cell based trends of water cover and spectral greenness. We further used structural equation modeling (SEM) taking into account trends of climate change and topographic conditions to infer how changes in surface water cover and vegetation could interact with each other to explain trends in spectral greening and browning. To account for the potential uncertainties associated with remotely sensed data, trend quantification methodology and spatial heterogeneity, we systematically tested for the robustness to different VIs, different trend indicators, different continents, different biome types, different permafrost conditions, and different observational scales.

Study area
We defined the study area using the biome map from World Wildlife Fund (Olson et al 2001). We extracted the polygons defined as 'tundra' and 'boreal forest/Taiga' in the WWF biome classification system as our study system (45 • N-80 • N). Before the analyses, we filtered out the areas that are subject to strong human activities using the Global Human Influence Index (HII) version 2 dataset at 1 km resolution (Wildlife Conservation Society-WCS and Center for International Earth Science Information Network-CIESIN-Columbia University 2005). The HII data values range from 0 to 64 representing an increasing intensity of human activities. We excluded the 1 km grid cells using a HII cutoff of 10, below which areas presenting apparent human land uses and areas in proximity to major roads and railways are mostly excluded, therefore human influences are considered very weak. We excluded the grid cells that intersect with coastlines and major rivers to avoid potential hydrological influences of these large water bodies. We also excluded the burnt areas during 2000-2015 using the MODIS MCD45 Burnt Area data (Giglio et al 2015) to avoid potential influences of fire.

Data sources and pre-processing
We used satellite-derived VIs to quantify vegetation dynamics at the circum-Arctic scale, including NDVI and EVI from the Terra Moderate Resolution Imaging Spectroradiometer (the MODIS MOD13Q1v006 level 3 product generated every 16 d) at a spatial resolution of 250 m (Didan 2015), and the Landsat images at a spatial resolution of 30 m. We adopted these multiple datasets to account for potential bias from single VIs or sensors, especially for highlatitude regions (Myers-Smith et al 2020). We selected the MODIS and Landsat data that are available for the study time period 2000-2015. For the MODIS NDVI and EVI data, we first picked up the pixels flagged as 'Good Data' in the Summary Quality layer, and then masked out the water pixels using the MODIS yearly Land Water Mask (MOD44Wv006) product (pixels classified as water in any year between 2000-2015). The Landsat NDVI and EVI were calculated from the atmospherically corrected USGS Landsat 7 Surface Reflectance Tier 1 dataset, where the pixels flagged as 'cloud' , 'cloud shadow' or 'snow' in the pixel_qa band were excluded. We also excluded the Landsat pixels classified as water in any year between 2000 and 2015 using the Global Surface Water Dataset (Pekel et al 2016). For both MODIS and Landsat VIs, we also masked out the non-vegetated pixels by excluding the pixels consistently presenting VI value below 0.1 across 2000-2015. Through these exclusions we are able to focus on quantifying spectral greenness of terrestrial areas. We then calculated annual mean growing-season (between 1 July and 1 September (Goetz et al 2005)), NDVI and EVI per vegetated pixel between 2000 and 2015 as indicators of spectral greenness. To test for the robustness with respect to observational scale, we additionally checked the 3rd generation NDVI data from the Global Inventory Monitoring and Modeling System (GIMMS NDVI3gv1, available for 1981-2013) at a spatial resolution of 8 km for the time period of 2000-2013.
The water cover data were obtained from the Global Surface Water Dataset (Pekel et al 2016) derived from Landsat data at 30 m spatial resolution. For the calculation of water cover of each year, we focused on the water bodies that can be consistently identified by all available cloud-free remotely-sensed data in the ice-free season (referred to as 'permanent water' in the Global Surface Water Dataset).

Spatial and statistical analyses
We divided the study area into 1 km × 1 km grid cells to calculate the mean yearly growing-season NDVI/EVI and percentage cover of permanent surface water within each grid cell. For the subsequent analyses we excluded the grid cells without any surface water cover throughout the study time period of 2000-2015, as we focus on the relationships between surface water cover and vegetation. We used a robust regression approach with the Theil-Sen estimator to determine the slope of the regression as an indicator of temporal trend (on an annual basis) for each variable (i.e. NDVI/EVI, surface water cover, MAT and MAP) per grid cell (figures 1(a), (b) and S1 (available online at stacks.iop.org/ERL/16/104004/mmedia)). We used the non-parametric Mann-Kendall trend test to test if there is significant monotonic increasing or decreasing trend. Considering that noise in the time series may potentially bias the trend estimates, we looked at two additional trend indicators, including signal-to-noise ratio (defined as trend divided by standard deviation of the time series), and significant trend (i.e. retaining the grid cells with p < 0.1 from the Mann-Kendall trend test for trends of VIs and surface water cover). We illustrated the distributions of mean NDVI/EVI trend within the state space of MAT/MAP trend against surface water cover trend (at bins sized of 0.003 • C yr −1 MAT in figure 2(a) and 0.2 mm yr −1 MAP in figure 2(b)). To test if the observed patterns are affected by grain size, we conducted multi-scale analyses by going through all abovementioned procedures on the basis of 2 km × 2 km, 4 km × 4 km, and 8 km × 8 km grid cells, respectively. We chose a static grid approach instead of a sliding window approach here.
We conducted a spatial neighborhood analysis to detect spatial associations between ecosystem changes in focal areas and neighboring areas. To quantify the effect of spatial autocorrelation (SA) processes in shaping the detected associations, we constructed SA models assuming (a) no SA, (b) SA generated purely by extrinsic processes (i.e. where the associations between spectral greenness and surface water cover are completely caused by certain spatially structured environmental factors), (c) SA generated purely by intrinsic processes (i.e. where changes of spectral greenness or surface water cover at focal sites can 'spill over' to influence the neighboring sites, causing the associations), and (d) SA generated by both extrinsic and intrinsic processes, respectively (see supplementary materials section S1 and Teng et al (2018) for details). If the models assuming existence of intrinsic processes (models 3 or 4) present better performance (assessed by the Akaike information criterion (AIC)), then there plausibly exist spatial interactions between spectral greenness and surface water cover changes. Therefore, these spatial modeling analyses can not only account for the potential influence of SA on model parameter estimation, but also help to infer the potential spatial interactions (Teng et al 2018).
We used piecewise SEM to explore potential drivers of the opposite trends of spectral greenness and surface water cover. In the SEM, we included MAT and MAP trends as climatic drivers of spectral greenness and surface water cover trends, that play a major role on how vegetation-water dynamics could respond to a long-term trend of global warming. Because topography affects drainage, we included topographic slope as an influencing factor on surface water cover, and therefore as an indirect driver of vegetation change. Such a minimal SEM is not meant to approach complex causal mechanisms in realistic situations. Instead, with this simplistic consideration, we aimed at testing if the remotely sensed observations are best fit under which one of the three different assumptions of interactions: (a) assuming one-way effect from spectral greenness trend on surface water cover trend, (b) assuming one-way effect from surface water cover trend on spectral greenness trend, and (c) assuming spectral greenness and surface water cover trends affect each other. As piecewise SEMs cannot disentangle feedback loops explicitly, bi-directional relationships are evaluated through simple correlation reflecting that the variables are co-varying. Corresponding with the abovementioned assumptions, we constructed three SEM models with (a) surface water cover trend as responsible variable, (b) spectral greenness trend as responsible variable, and (c) covarying surface water cover trend and spectral greenness trend. We used simple linear models in the piecewise SEM as the data approach normal distributions. We assessed AIC based on the Fisher's C statistic as the indicator of model performance. The models were fitted for the whole tundra and boreal biomes, as well as for different continents, different biome types, different permafrost conditions, different VIs and different trend indicators separately. Data processing and spatial analyses were conducted with the Google Earth Engine platform (Gorelick et al 2017) and Arc-GIS 10.2. The SEM was conducted using the R 3.5.3 software (R Core Team 2018) with the piecewiseSEM package (Lefcheck 2016).

Divergent trends of surface water cover and vegetation
Our analyses show that spectral greening and browning co-occur across high latitude regions where water bodies are present (figure 1). Both spectral greening and browning can also happen as MAT (figure 2(a)) and precipitation change ( figure 2(b)). These puzzling patterns reflect the complexity of local tundra and boreal vegetation dynamics in response to climate changes (Myers-Smith et al 2020).
Our analysis demonstrates an interesting spatial congruence between opposite trends in surface water cover and vegetation condition reflected by different VIs analyzed (figure 1(c), table S1). This negative correlation becomes even stronger if areas with insignificant temporal changes in these variables are excluded (Mann-Kendall trend test, p > 0.1), as seen from a drastic increase of Pearson's correlation coefficients from ∼0.2 to ∼0.6 (table S1).
By slightly smoothing the VI trends, a clear pattern emerges, namely, we find spectral greening to be strongly associated with reductions in surface water cover, whereas spectral browning is strongly associated with increases in surface water cover (figures 2(a) and (b)). These patterns of greening and browning in relation to opposite trends in surface water cover are highly consistent for both NDVI and EVI derived from MODIS (250 m resolution) and Landsat (30 m resolution) data, and are also robust across regions with contrasting permafrost conditions (i.e. continuous permafrost, discontinuous permafrost, and permafrost-free areas), biome types (i.e. boreal and tundra biomes), and continents (i.e. Eurasia and North America) (figure S2). Interestingly, the patterns are most pronounced at smaller grid cell sizes of 1 and 2 km, and become attenuated at 4 and 8 km ( figure S3). Two landscape-scale examples of these patterns can be seen at the Kytalyk site in eastern Siberia, where ecosystem changes are characterized by spectral browning associated to surface water cover increases, and in the Atqasuk site in North America, where spectral greening is associated to decreases in surface water cover (figure S4).
These revealed associated opposite trends between surface water cover and spectral greening and browning could partly be a result of reciprocal (a) Map of spectral greening and browning reflected as increasing and decreasing Landsat EVI trends respectively at 1 km grid cells. Areas subject to human influences, fire and low-quality of remote sensing data were filtered out (see section 2). (b) Map of increasing and decreasing surface water cover trends calculated using the Global Surface Water Dataset (Pekel et al 2016) derived from Landsat data. (c) Spatial congruence between spectral greenness and water cover trends (Pearson's r correlation coefficient of −0.18, p < 2.2 × 10 −16 , n = 5191 139). The point density and the slope line based on ordinary least square regression are shown (note that the 95% confidence interval is very close to the black solid slope line). The correlation becomes stronger (Pearson's r correlation coefficient of −0.60, p < 2.2 × 10 −16 , n = 6886) if looking at areas with significant trends only. See figure S1 and table S1 for the results based on all four VIs used (i.e. Landsat NDVI and EVI; MODIS NDVI and EVI).

Figure 2.
Spectral greenness trends in relation to surface water cover and climatic trends. (a) Distribution of spectral greenness trend in relation to MAT trend as x axis and surface water cover trend (as y axis). (b) Distribution of spectral greenness trend in relation to MAP trend as x axis and surface water cover trend as y axis. (c) A minimal SEM suggests plausible interactions between the trends of surface water cover, spectral greenness and climate. Both spectral greening (positive EVI trend represented by blue dots) and browning (negative EVI trend represented by red dots) can happen when the climate has become warmer or cooler, drier or wetter ((a) and (b)). Interestingly, spectral greening mostly occurs in areas with decreasing surface water cover, whereas browning mostly occurs in areas with increasing water cover ((a) and (b)). The horizontal dashed lines in ((a) and (b)) separate between decreasing surface water cover (negative water cover trend) and increasing surface water cover (positive water cover trend). The model assuming co-varying surface water cover and spectral greenness trends (a bi-directional relationship) is the best fit across the tundra and boreal biomes on a global scale (c). Red and blue arrows indicate negative and positive effects (with their standardized coefficients, p < 2.2 × 10 −16 for all effects), respectively. Dashed double-direction red arrow indicates co-variation of spectral greenness trend and surface water cover trend, possibly resulting from their reciprocal feedbacks. For model fitting, bootstrapping was conducted with 1000 repetitions and 50 000 random samples drawn from 51 91 139 data points in each repetition to account for potential influence of SA on parameter estimates. Note here that p values greater than 0.05 indicate satisfactory model performance. Empirical evidence has also been obtained from localscale experiments in northeastern Siberian tundra showing that where shrub canopy has been removed ponds quickly form (Nauta et al 2015).
Interestingly, our further analysis reveals that when a particular focal area experienced a decreasing trend of surface water cover, its neighboring terrestrial areas tended to exhibit spectral greening (figures S5-S7). The association of changes was also evident in the opposite direction as spectral browning increased adjacent to newly formed water bodies (figures S5-S7). These neighborhood associations suggest the existence of 'spillover effects' (i.e. a spatial process by which changes in a particular area can cascade to induce changes in nearby areas) (LeSage and Pace 2009). We used spatial regression modeling to infer if such neighboring associations could be purely attributed to a SA process (i.e. where certain external environmental factors drive both greenness and surface water cover trends, thus giving rise to their associations, a process referred to as extrinsic SA (Teng et al 2018); see section 2 and supplementary materials section S1). Our results suggest that the observed associations are not solely caused by extrinsic SA. Instead, we found signs that changes in both spectral greenness and surface water cover at focal sites can spill over to produce significant influences on land cover changes in neighboring areas (figure S8), suggesting that vegetation condition and surface water cover plausibly interact with each other.

Inference of vegetation-water interactions
The SEM suggests that overall the associated trends of spectral greenness and surface water cover are better explained by bi-directional relationships, in which spectral greenness and surface water cover trends influence each other. In some situations, the performance of this bi-directional model may not exceed the results assuming one-way effects of surface water cover on spectral greenness (see figure 2(c) for the model structure and table S2 for results depending on continent, vegetation type and permafrost condition, as well as type of vegetation index and trend indicator used). We also observed a negative effect of topographic slope on water cover trend (figure 2(c)), in line with the well-known effect of slower drainage in flat areas facilitating water accumulation and increases in surface water cover. However, the effects of topographic condition and climate (i.e. MAT and MAP) trends, on changes in surface water cover, were much weaker than the inferred vegetation-water interactions, echoing the highlighted complexity of soil and surface water changes in response to climate warming and permafrost thawing in northern high latitudes (Walvoord and Kurylyk 2016).

Possible mechanisms
These revealed associations between changes in spectral greenness and surface water cover could be explained by several mechanisms. For example, increasing waterlogged conditions resulting in increases in surface water cover can elevate mortality of trees and shrubs as reported in lowland boreal forests and tundra across Alaska, Canada, and northern Eurasia (Forbes et al 2010). Loss of woody cover may exacerbate summer thawing of permafrost, further increasing surface water cover (Blok et al 2010, Nauta et al 2015. Interestingly, bi-directional relationships appear to be mostly associated to regions with continuous permafrost where soil ice content is high (table  S2). This is in line with the suggested feedback between permafrost thawing and vegetation cover loss (Blok et al 2010, Nauta et al 2015. If permafrost condition indeed underlies the direction of vegetationwater interactions, one would expect to observe this bi-directional pattern more commonly in regions where permafrost is more abundant. In line with this, we found more pronounced bi-directional relationships in the tundra biome dominated by continuous permafrost than in the more southern boreal biome (table S2). We also found more pronounced bi-directional relationships between trends in surface water cover and spectral greenness in Eurasia (whereas one-way in North America, as suggested by the SEM fitting). This continental difference could be attributed to the larger fraction of continuous permafrost in Eurasia than in North America (Brown et al 2002) as well as the stronger heatwaves occurring in Eurasia (Schubert et al 2014).
Despite the dominant negative correlation between trends of spectral greenness and surface water cover, it is not uncommon to find positive associations in particular areas scattered across the boreal and tundra biomes. In relatively dry areas within a landscape, increasing surface water can alleviate limitation of soil water and facilitate vegetation growth (Ruiz-Pérez andVico 2020, Sungmin et al 2020). Also, in wetlands, the productivity of wetlandadapted plants may increase with warming and CO 2 fertilization, resulting in greenness increases along with increases in surface water cover induced by permafrost thawing (Park et al 2016). On the other hand, drier soils resulting from surface water drainage (Smith et al 2005, Hinkel et al 2007, Marsh et al 2009 can facilitate the replacement of sedges and mosses by shrubs and trees (Frohn et al 2005, Jorgenson et al 2013, Li et al 2017 and enhance tree growth (Ropars et al 2015). In turn, positive interactions between shrubs and trees can reinforce woody plant encroachment and succession , Limpens et al 2021, increasing water transport from soil to atmosphere and further contributing to reductions in soil water and surface water cover (Waddington et al 2015).
While soil moisture has been shown to correlate with Arctic greening (Berner et al 2020), the explicit link between surface water and vegetation change has been poorly explored using remote sensing as surface water bodies are often excluded to prevent confounding factors. Given the significant complexity and heterogeneity of ecosystem processes in northern high latitudes, we do not expect that the inferred interactions between vegetation and surface water changes are ubiquitous at the biome-wide scale. However, they may at least partly explain why both greening and browning have been observed across the boreal and tundra biomes. These interactions between vegetation condition and water cover may alter carbon dynamics and biogeochemical cycles. Vegetation browning may result in loss of stored above-ground and below-ground carbon as permafrost thaws and stored organic matter decomposes (Walter et al 2006, Schuur et al 2015. In well drained areas, reduction of surface and soil water may also lead to enhanced decomposition and increase CO 2 and methane release (Natali et al 2015). On the other hand, while vegetation greening can partly compensate for the loss of stored carbon, it may also further stimulate decomposition of native soil carbon (Fontaine et al 2003).

Caveats and limitations
We cannot rule out the possibility that intrinsic limitations of current remote sensing tools may influence the observed patterns to some extent. For example, 'gridding artifacts' (Tan et al 2006) denote that nearby water pixels can contribute to the reflectance of the focal vegetated pixel. Sufficiently strong gridding artifacts thus can result in spectral browning positively correlated with increasing cover of surface water. However, if the gridding artifacts were fully responsible for the surface water-browning pattern, one would expect to observe clear differences between MODIS and Landsat data. Yet, both remote sensing tools detect the same patterns (figure S3).
Another technical caveat to keep in mind is the influence of subpixel heterogeneity, which is inevitable at medium to coarse spatial resolutions. The probability that focal pixels are partly covered by surface water probably increases and leads to spectral browning if neighboring pixels exhibit a positive trend of surface water cover and/or are largely covered by surface water. However, if the observed negative trend associations were largely an artefact resulting from subpixel heterogeneity, Landsat data that are less prone to subpixel heterogeneity would have given weaker trend correlations than the MODIS data. Yet, the results from MODIS and Landsat data are consistent with each other (table S1).
Despite that our observed patterns are unlikely resulting from pure remote sensing artifacts, the abovementioned caveats should be taken into account for further unraveling vegetation-water interactions in the northern high latitude regions using remote sensing data. The increasing availability of high-resolution remote sensing datasets at multiple temporal coverages and enhanced computation power of cloud computation platforms may help overcoming the current technical limitations of coarse-and moderate-resolution remote sensing tools.

Conclusions
Taken together, our analyses reveal clear associations between divergent changes in surface water cover and spectral greening and browning. These changes we infer may also create positive feedback loops that selfreinforce changes in vegetation, permafrost thawing and surface water at relatively small spatial scales. In local areas where these feedbacks are strong enough, they may effectively reinforce or offset the direct effects of climate changes on permafrost thawing and vegetation. While it is impossible to rigorously verify the existence and quantify the intensity of these feedbacks globally with correlational evidence only, our results convey the important message that changes in surface water could play a critical (but largely neglected) role in ecosystem structure and functioning at northern high latitudes. Our results thus imply that the recent projections of abrupt permafrost degradation across northern high latitudes (Teufel and Sushama 2019) can be locally reinforced by local interactions between water and vegetation (Scheffer et al 2012). This may help to explain the existence of abrupt vegetation states reported earlier across the boreal and tundra biomes (Scheffer et al 2012, Xu et al 2015, Teufel and Sushama 2019. Searching for abrupt changes in VI time series using methods such as the breaks for additive seasonal and trend method (Verbesselt et al 2010) may help to better foresee if and to what extent biome-wide tipping points could occur in the future.
Our findings have major implications globally. Climate warming is expected to increase permafrost thawing. The accumulation of inland surface water as a consequence, especially in Arctic and sub-Arctic regions with low drainage, will likely increase the extent of wetlands, changing ecological communities, biogeochemical cycles and their feedback to the climate system. This process could be accelerated in regions where warming is also combined with higher rainfall levels or stronger pulses of precipitation events as our planet warms up. We expect therefore that vegetation browning may be amplified in many regions across the northern high latitudes. To monitor the changes we expect and assess the mechanisms we propose here, it will be important to combine high resolution remote sensing with field experiments at larger spatial extents that currently undertaken to unraveling the complex interactions between climate, vegetation cover, snow accumulation, permafrost condition, hydrology and albedo. Simulation models that adequately account for these interactions may help to better project the responses of these massive ecosystems to climate change.

Data availability
All data used in this study are publicly available. The MODIS NDVI and EVI data are available at https://lpdaac.usgs.gov/products/mod13q1v006/. The Landsat NDVI and EVI data are available at https://developers.google.com/earth-engine/data sets/catalog/LANDSAT_LE07_C01_T1_SR. The water cover data are available at https://globalsurface-water.appspot.com/download. The MAT and MAP data are available at the CRU (www.cru.uea.acuk/data). The permafrost data are available at https://nsidc.org/data/ggd318. The elevation data are available at www.usgs.gov/centers/eros/ science/usgs-eros-archive-digital-elevation-global-30-arc-second-elevation-gtopo30.
All data that support the findings of this study are included within the article (and any supplementary files).