Climatic Drivers of Greening Trends in the Alps

: Since the 1980s, vegetated lands have experienced widespread greening at the global scale. Numerous studies have focused on spatial patterns and mechanisms of this phenomenon, especially in the Arctic and sub-Arctic regions. Greening trends in the European Alps have received less attention, although this region has experienced strong climate and land-use changes during recent decades. We studied the rates and spatial patterns of greening in an inner-alpine region of the Western Alps. We used MODIS-derived normalized difference vegetation index (NDVI) at 8-day temporal and 250 m spatial resolution, for the period 2000–2018, and removed areas with disturbances in order to consider the trends of undisturbed vegetation. The objectives of this study were to (i) quantify trends of greening in a representative area of the Western Alps; and (ii) examine mechanisms and causes of spatial patterns of greening across different plant types. We show that 63% of vegetated areas experienced signiﬁcant trends during the 2000–2018 period, of which only 8% were negative. We identify (i) a climatic control on spring and autumn phenology with contrasting effects depending on plant type and elevation, and (ii) land-use change dynamics, such as shrub encroachment on abandoned pastures and colonization of new surfaces at high elevation. Below 1500 m, warming temperatures promote incremental greening in the transition from spring to summer, but not in fall, suggesting either photoperiod or water limitation. In the alpine and sub-alpine belts ( > 1800 m asl), snow prevents vegetation development until late spring, despite favorable temperatures. Instead, at high elevation greening acts both in summer and autumn. However, photoperiod limitation likely prevents forested ecosystems from fully exploiting warmer autumn conditions. We furthermore illustrate two emblematic cases of prominent greening: recent colonization of previously glaciated/non vegetated areas, as well as shrub/tree encroachment due to the abandonment of agricultural practices. Our results demonstrate the interplay of climate and land-use change in controlling greening dynamics in the Western Alps.


Introduction
Since the 1980s, vegetated lands have experienced widespread greening at a global level [1]. Growing season integrated leaf area index has increased significantly across 25-50% of the vegetated areas [2], and the normalized difference vegetation index (NDVI) has increased globally at a rate of 0.0013 yr −1 for the 1982-2011 period [1]. Accurate descriptions of greening dynamics and controls are crucial because increased greening is associated with increased carbon uptake by vegetation [3], with important implications on the carbon dioxide balance [4]. A large body of literature shows that Arctic and sub-Arctic environments are greening faster than temperate areas in response to accelerated warming in these regions [5][6][7][8]. However, substantially less work has been devoted to quantifying and understanding patterns of greening in temperate mountain regions, such as the European Alps. Yet, the Alps represent a key area where climate and land-use changes are already strongly influencing ecosystem processes, and shifts in vegetation structure have immediate consequences for services provided to densely populated alpine valleys. Temperature increase is about double than the global average [9,10], and snow cover duration has decreased above 2000 m in the last 50 years [11], although climatic trends since the 1990s are weaker [12]. New surfaces are becoming available for primary succession, including recently deglaciated areas [13], and an upward shift of alpine vegetation has been observed, documenting the colonization of previously non-vegetated areas and increasing species richness on mountain summits [14]. Other surfaces are radically changing because of a combination of land-use and climatic change, resulting, for example, in tree-line rise [15,16], shift in species [14,17], and pasture abandonment. A recent study [18] documented widespread greening in grassland, low-shrub, and scree slopes of a protected area in the French Alps and proposed rising temperatures and decreased snow-cover duration as the main drivers. However, the positive effect of rising temperatures can be modulated by other factors (e.g., drought, photoperiod) in a complex inter-play. For example, water deficit may offset the benefit derived from warmer temperatures in mountain grasslands [19]. Forested ecosystems could also be subjected to greening trends via biomass increase, forest densification, and shifts in species composition. An alpine-scale study on land surface phenology showed extensive lengthening of the growing season and demonstrated the role of snow duration on phenology trends in the last 20 years [20].
Changes in vegetation structure can be attributed to two categories: (1) gradual changes that may determine increase or decrease of plant productivity, generally referred to as greening or browning [21], and (2) abrupt changes due to disturbances, such as fire, windstorms, avalanches, or land reshaping. Indeed, de Jong et al. [21] demonstrate the importance of accounting for abrupt and gradual changes when analyzing long-term NDVI time series. To the best of our knowledge, abrupt changes have not been taken into account while studying greening trends in the Alps. Furthermore, Emmett et al. [22] showed the importance of disentangling climate and disturbance effects when studying greening trends. One of the most prominent tools for the analysis of trends and abrupt changes in time series of vegetation indices is the Breaks For Additive Seasonal and Trend (BFAST) method [23,24], which was successfully used in global and regional analyses. In particular, a recent study [25] deployed the BFAST method on MODIS-derived NDVI time series in a mountain environment and confirmed the robustness of this approach also in complex terrain.
In this study, we used 19 years of MODIS-derived NDVI to answer the following questions: (1) What is the magnitude of greening trends in the Western Alps when the role of the disturbance is accounted for? (2) What are the mechanisms and causes of the spatial patterns of greening across plant types, elevation gradients, and seasons?
We hypothesize that (1) greening trends at high elevation documented for other Alpine regions may be due to warmer temperatures and longer snow-free periods; and (2) warmer summer temperatures in the valley lowlands can also lead to greening, and the BFAST method can help to identify it.

Study Area, Climate Data, and Vegetation Mapping
This study was conducted in Aosta Valley, an inner-alpine valley located in north western Italy. The large elevation range between 300 and 4800 m (mean: 2100 m), over a relatively small surface (3262 km 2 ), produces steep temperature gradients, with mean annual temperatures of 14 • C at 500 m down to −3 • C at 3000 m for the period 2000-2018. The continental character of the valley results in pronounced precipitation gradients, with cumulative values of more than 1800 mm year −1 in peripheral areas (west and south-east boundaries), and local minima of 500 mm year −1 in the center of the region.
Daily gridded temperature, precipitation, and evapotranspiration data at 100 m spatial resolution were used as climatic descriptors for the period 2000-2018. Precipitation and temperature are derived from geostatistical interpolation of ground data from the regional network of automatic weather stations [26]. Daily temperature (precipitation, evapotranspiration) maps were averaged (cumulated) to a long-term annual product, as well as seasonally, considering meteorological seasons. A simplified water balance was estimated as the difference between precipitation and evapotranspiration. Evapotranspiration was computed using the hydrological model Continuum [26,27]. The model is distributed and solves the conservation equation for mass and energy at the soil on a regular grid with 100 m resolution. It estimates latent and sensible heat fluxes using a bulk transfer approach based on bulk transfer formula [28]. It is used operationally in the region for flood forecasting and water managements. Maps of snow water equivalent (SWE, 500 m spatial, 15 days temporal resolution) were obtained from a statistical model combining MODIS snow covered area, snow depth data from a ground meteorological network, manual SWE measurements, and topographical attributes in a multiple-regression-krieging framework. These data were used to describe the snow climatology of the study area, by means of two average maps: mean maximum SWE and mean snow melt date for the period 2002-2018.
Vegetation types were mapped by photo-interpretation of aerial photographs collected in the years 2005 and 2006. Vegetation mapping was performed according to the habitats defined in Annex I of the Habitats Directive [29] and described in detail in [30]. The final product, digitized in the form of polygons at a scale ranging between 1:5000 and 1:10,000, was then rasterized at 250 m to match the spatial resolution of the satellite products. The 56 habitats identified were further aggregated in 8 plant functional types (PFT, Figure 1). Forests cover 34% of the study area, with a 20% contribution of deciduous species and 14% of evergreen needleleaf forest. The most abundant species is larch (the only deciduous needleleaf species, DNF), which alone covers 14% of the whole surface. Grasslands cover 26% of the surface, with only 4% meadows and 21% pastures. The second most abundant type is vegetated scree slopes, with 29% cover. Finally, 5% of the surface is covered by glaciers and recently deglaciated areas.

Satellite Data
All satellite data cover the period from February 2000 to December 2018. MODIS Terra surface reflectance images (MOD09Q1) collection 6 at 250 spatial and 8 days temporal resolution [31] were downloaded from National Aeronautics and Space Administration (NASA) ftps using the MODIStsp R package [32]. The same was used to compute NDVI maps, which were then masked according to the MOD09Q1 quality flags; only pixels classified as of "best quality" were retained. NDVI was preferred over the enhanced vegetation index (EVI) for consistency with the majority of studies on global greening trends, which considered NDVI only. The MOD10A2 collection 6 snow covered area (SCA) standard product was also downloaded and re-sampled at 250 m spatial resolution using the nearest neighbor interpolation method. The SCA dataset was used to compute a winter baseline for NDVI with a method similar to [33]. First, the resampled SCA dataset was used to mask snow-free pixels. For each pixel, the 5th percentile of NDVI was computed on the whole time series to determine a per-pixel climatological minimum snow-free NDVI. This minimum snow-free map was used to mask each 8-day NDVI map, with values lower than the climatological minumum substituted with it on a pixel basis.
Corrected 8-day maps were averaged to a monthly product which was used in all subsequent analyses.

Data Analysis
The BFAST (Breaks For Additive Seasonal and Trend) method implemented in the bfast R package was used to calculate NDVI trends over the 2000-2018 period (TREND BFAST ). bfast integrates the decomposition of time series into trend, seasonal, and remainder components with methods for detecting and characterizing abrupt changes within the trend and seasonal components [23,24]. As such, this method is particularly suitable for the detection of disturbed time series and concomitant trend detection. Pixels showing clear breaks in either long-term trend or seasonal components were removed from the analysis of trends.
Because the wide range in mean NDVI across different PFTs can complicate the interpretation of absolute NDVI trends, TREND BFAST values were also expressed as percentage of the long term pixel-based mean and further multiplied by 20 to represent the cumulated NDVI increase in the ∼two MODIS decades 2000-2018 (see Figure 2b). We also calculated another metric derived from NDVI, the area under the the curve (AUC, i.e., the NDVI integral), used as a surrogate for aboveground productivity [34]. AUC was computed after subtracting the minimum baseline to the seasonal values. In addition to the BFAST method, we computed NDVI trends with the approach proposed by [18]. We extracted maximum summer NDVI for each year and grid cell and calculated the trend by the Theil-Sen median slope method and Mann-Kendall method for significance (TREND MAX , [35]). Pixel-based trends for temperature and precipitation data were computed by decomposing the time series according to [36]. The slope of the trend component was then estimated with the Theil-Sen median slope method [35]. Precipitation trend was rarely significant and discarded from subsequent analysis.
The relationship between TREND BFAST and environmental/climatic variables was investigated with Random Forests (RF) using the ranger package [37]. Explanatory variables included in the model are summarized in Table 1. Several important climatic variables were excluded because of high autocorrelation, identified by inspecting variable importance table: these include average temperatures (e.g., annual T, spring T both correlated to summer T) and phenological variables (growing season length, snow melt date, both correlated to summer T and maximum SWE). All models were run with the caret package [38] in R [39]. Model performance metrics are calculated from 20-fold cross-validation. As for RF parameters, we let vary the number of variables to possibly split at in each node (mtry) and the splitting rule (variance vs. extratrees, see [40]), whereas minimal node size and tree number were kept as default (5 and 500, respectively).
The dynamics of colonization of new surfaces was showcased by extracting NDVI time series in the proglacial area of the Rutor Glacier (45 • 40 12 N, 6 • 59 11 E, 2560 m a.s.l.). We focused on three points corresponding to markedly different ages of deglaciation, which was defined by visual inspection of historical aerial photographs retrieved from the web map service (WMS) of the Geoportale Nazionale (http://www.pcn.minambiente.it/mattm/en/view-service-wms/). Table 1. Environmental and climatic variables used in the random forest (RF) model.

Variable Type Details Acronyms
Thermal conditions Mean summer temperature, temperature trend, heat load index [41] Tsummer, Ttrend, hli Water status Summer water balance, maximum snow water equivalent, annual and summer precipitation PmET, maxSwe (snow water equivalent), P, Psummer Vegetation type plant macrotype as in Figure 1 and mean max normalized difference vegetation index (NDVI) pft, maxNDVI All analyses were conducted with the R software for statistical computing, version 3.6.1 [39].

Comparison of TREND BFAST and TREND MAX
The BFAST method applied to the 2000-2018 NDVI imagery on Aosta Valley revealed that 31% of pixels showed a disturbed time series characterized by one or more break points (Table 2). Disturbed pixels were mainly located at low elevation (69% below 1500 m). Deciduous forests showed the highest percentage of disturbed pixels, but evergreens and meadows followed closely. Because an analysis of the break points was beyond the scope of this paper, disturbed pixels were masked out before computing the final trend map (Figure 2). A small fraction (6%) of the non-disturbed pixels showed no trends, while the remaining (63%) were characterized by a significant NDVI trend. Among these, only 8% showed negative trends, whereas 92% of trends were found to be positive. TREND BFAST ranged between −0.001 yr −1 and 0.004 yr −1 and was slightly positively skewed (mean = 0.0011, median 0.00094). Different plant functional types were characterized by distinct absolute trend values; deciduous broadleaf forests (DBF) and meadows had the highest rates of greening, followed by evergreen and larch forests (DNF). Relative trends (expressed as percentage of the mean NDVI and integrated over the length of the time series, Figure 2c,d) ranged between ∼−9 and 55%. The highest median values above 30% were found for Glacier Forefield category. Scree slopes vegetation also showed a median increase of greening higher than 15%. For the other plant types, relative TREND BFAST were rather homogeneous. TREND MAX were on average three times higher than correspondent TREND BFAST (Table 3). However, the number of pixels with significant TREND MAX were only 16%, much lower than for TREND BFAST (63%). In addition, we found a complex relationship between trends retrieved with the two methods ( Figure 3). The expected linear relationship between the two trends holds as long as maximum summer NDVI remains below 0.6. At higher summer NDVI values, only the BFAST method was able to further detect increasing trends, and the relationship between the two trends reversed. Correspondingly, the number of significant pixels detected by the TREND MAX method decreased. This behavior suggests that the TREND MAX method is more prone to NDVI saturation compared to TREND BFAST method. Table 2. Proportion (%) of pixels showing a significant trend or disturbed time series based on the Breaks For Additive Seasonal and Trend (BFAST) method, separated according to plant macro-types and two elevation bands. The proportion is calculated by dividing the number of disturbed pixels by the total number of pixels of a given plant type, except for the totals in parentheses (bottom row), which refer to the number of pixels in the elevation band.

Factors Controlling NDVI Trends
The RF model showed good performance in describing NDVI trends (TREND BFAST ) as a function of environmental/climatic drivers (R 2 (CV) = 0.50, RMSE (REL) = 35%, MAE (REL) = 26%). The relative importance analysis revealed that summer temperature was the most influential variable, followed by maximum snow water equivalent, and to a lesser extent, annual cumulative precipitation. The bivariate relationships between NDVI trend and the two most important climatic/environmental variables are shown in Figure 4. TREND BFAST is lower at colder locations and increases linearly with temperature, and it is more sensitive to temperature in the lower and upper temperature range, rather than in the median (6-12 • C) range. TREND BFAST is also inversely correlated to SWE. Summer temperature and maximum SWE are strongly correlated to each other (r = −0.9) and to elevation (−0.97 and 0.92 for temperature and SWE, respectively), showing that the interplay of thermal and snow control on vegetation greenness co-varies with elevation, and the latter can be used as a proxy that integrates both snow and temperature. To consider the evolution of NDVI in the period 2000-2018 under markedly different thermal and snow regimes, we analyzed separately pixels above 1800 m and below 1500 m. These elevation thresholds were chosen because on average all pixels below 1500 m are snow free on March 1st and we can assume no snow limitation on spring phenology. Pixels above 1800 m are snow-covered at list until May 1st, and we can hypothesize significant snow limitation on spring vegetation development above this elevation.
On these two subsets, we compared monthly NDVI trajectories for the two periods 2000-2004 and 2014-2018 as representative of the beginning and end of the study period (Figure 5i,j). NDVI seasonal profile has evolved very differently in the two elevation bands from 2000 to 2018. Above 1800 m, former and current NDVI seasonal patterns were very similar during the first part of the growing season, whereas current NDVI seasonal pattern is characterized by a higher maximum and a sustained greening throughout summer and early fall. Below 1500 m, differences between former and current NDVI trajectories are already evident in spring, and the difference tends to reduce in late summer and fall. The general pattern described above is valid for all PFTs with some peculiarities (Figure 5a-h panels). Forested ecosystems above 1800 m show weaker greening in autumn compared to grasslands, shrublands, and screes. The difference between former and current areas under the curve (AUC) is positive at all elevations and for all PFTs but differs in magnitude, with twice as much of an increase above 1800 m compared to below 1500 m. Scree slopes display the highest increase in AUC (40%) during the study period, followed by shrublands, grasslands, and evergreen forests. AUC in deciduous forests increased by an equal rate of 8.4% at all elevations.

Methodological Considerations
We studied the rates, spatial patterns, and potential controls of greening in an inner-alpine region of the Western Alps by analyzing MODIS imagery for the period 2000-2018. The BFAST method allowed us to remove disturbed and non-linear signals from ∼37% of the vegetated land and to focus our analysis on widespread greening trends (Figure 2) representing more than half of the vegetated areas. Compared to the commonly used TREND MAX method [18], TREND BFAST allows the exploitation of the whole time series [23] and thereby accounts for changes to seasonal NDVI profile that the TREND MAX method neglects. We demonstrated that this method marks an improvement over the traditional TREND MAX ; not only does it allow the detection of a greater number of pixels with significant trends (i.e., more statistical power) but it is also more robust against the saturation of NDVI signal (Figure 3). This issue, for example, prevented Carlson et al. [18] from investigating trends over forested areas with maximum NDVI higher than 0.75 in an Alpine region similar to the one analyzed in this study. Removing disturbed time series allowed us to exclude recent land-use change from this analysis; however, changes that occurred prior to year 2000 may still be affecting observed trends by means of legacy effects of past land-use practices or consistent land-use practices during recent years (i.e., extensive grazing). Finally, the moderate spatial resolution of MODIS data (250 m) limits our ability to interpret NDVI trends at small-scale, as, for example, at the level of alpine plant communities or in the case of changes in land-use or practices on limited surfaces.

Controls and Mechanisms of Greening
Snow and temperature regimes control the shape of greening in this region ( Figure 5). The recent increase in air temperature recorded all over the Alps results generally in an earlier start of the growing season [20]. In our study, we confirmed this pattern below 1500 m (Figure 5b), where snow melts before March 1st and plant development can take advantage of warmer spring temperatures. Warmer springs therefore result in sustained greening in the transition from spring to summer, whereas the same is not true in the transition to fall, suggesting either photoperiod [42][43][44] or water limitation [45]. Between the two, we hypothesize water could be more likely the limiting factor, because (i) if photoperiod were the cause, a differential response between forests (more sensitive) and grasslands would be expected, but this was not the case ( Figure 5), and/or (ii) lowlands in this inner alpine region receive rather low summer precipitation amounts (<600 mm/yr), and water availability could be critical while approaching the end of the growing season. This is confirmed by the cumulated summer/fall water balance (precipitation-actual evapotranspiration) calculated for the two elevation belts, that averages 346 ± 71 and 179 ± 62 mm above 1800 m and below 1500 m, respectively (mean ± 2SE for the period 2008-2017). Even though an average value of 179 mm below 1500 m cannot be considered indicative of water stress, it may still represent sub-optimal conditions, compared to those observed above 1800 m.
Above 1800 m the snow cover prevents the vegetation development until late spring. While a trend toward an earlier snow melt was reported at an alpine-wide scale by [20] for the 2000-2018 period, longer records suggest little evidence of earlier snow melt in the Alps after the 1990s [46]. For this reason, we hypothesize plants above 1800 m cannot take advantage of warmer springs due to the presence of snow cover [34]. It is likely, however, that this hypothesis applies only to the last 20 years, while before the 1990s, the combination of warmer temperatures and earlier snow melts determined an earlier beginning of the growing season even at high elevation, as documented by previous studies [47,48]. Instead, the majority of greening at high elevation is manifested toward summer and autumn, with remarkably higher NDVI values until October. Autumn plant senescence at high altitude occurs normally well before the snow onset, which indeed is not a limiting factor and the vegetation responds to warmer fall temperatures. An important role of photoperiod is suggested by the divergent behavior of trees versus grass-dominated ecosystems at high elevation (Figure 5a-d).
The late-season greening trends strikingly differed between high elevation non-forested vegetation and forests. The pronounced trend for grasslands may reflect a predominant control by temperature. In contrast, the weaker response of forest to late-season warming may reflect a predominant control by photoperiod [49,50].
Overall, the negative relationship between greening trend and maximum winter SWE (representing the quantity of water stored in the form of snow available in the watersheds at snow melt, Figure 4) seems counter-intuitive because we would expect that more water available may result in higher greening rates. However, we hypothesize that the observed inverse relationship reflects the control that snow cover exerts on spring phenology through the decoupling between air, canopy, and soil temperatures. Higher maximum SWE values are in fact associated with later snow melt, and both these conditions are more likely to occur at higher elevation, where trends are lower in magnitude. These results suggest that the water input from snow melt plays a marginal role in regulating the phenology of alpine ecosystems.

Magnitude of Trends and Spatial Patterns
Large differences in relative greening across plant types suggest that different mechanisms may be responsible for observed trends. NDVI in glacier forefields increased by up to ∼50% over the 2000-2018 period because of the colonization of previously glaciated areas [51,52]. To illustrate this feature, we selected MODIS pixels located in the proglacial area of the Rutor Glacier, one of the largest glacier bodies in Aosta Valley (Figure 6), based on historical ortophotographs. Areas deglaciated well before 1988 (at least 10 years) show a weakening NDVI trend, whereas the fastest greening rates are apparent at location where primary succession started ∼30 years ago. Recently deglaciated areas (i.e., less than 20 years ago) are greening at an intermediate rate. This qualitative analysis shows that MODIS spatial resolution may be sufficient to inspect the phenomenon of primary succession in areas with flat and homogeneous proglacial morpohology. The process could be further analyzed using data from different satellites, like LANDSAT, providing longer time series at higher spatial resolution. In line with findings by [18], scree slopes at high elevation show a greening increase of more than 10% over the 2000-2018 period. This increase can be associated to plant colonization of previously non-vegetated surfaces consistent with the observed upward shift of herbaceous plant species [14,53]. Shrub and tree encroachment may be responsible for the ∼10% greening observed on shrublands, grasslands, and meadows. This process results largely from the abandonment of agricultural practices in previously managed areas [54] and can be in part attributed to rising temperatures, as well [55]. The magnitude of this process is further testified by the AUC increase of grasslands at high elevation, which is more than double than at low elevation, due to the fact that shrub encroachment caused by pasture and agricultural practices abandonment occurs more often at higher elevation than in the lowlands [56]. These phenomena are often localized and difficult to detect at 250 m spatial resolution. However, we illustrate the process in Figure 7. The visual comparison documents a sharp change in vegetation, reflected by a significant trend (0.004/yr, 15% /20 yr). Maximum summer NDVI in this area increased from 0.66 in 2000 to 0.85 in 2018. Figure 7. Shrub/tree encroachment illustrated by successive aerial photographs on a former meadow. TREND BFAST is annotated. MODIS grid is drawn for reference. The pixel for which values are annotated is drawn in red.

Conclusions
This study provides evidence for widespread greening in the Alps and explores the climatic cues controlling its magnitude and spatial patterns. We found that warmer temperatures promote an increase in plant productivity, but that the seasonal timing of greening is regulated by other elevationand ecosystem-dependent factors. Among these, snow prevents spring plant development above 1800 m and photoperiod likely limits forested ecosystems in autumn, whereas at low elevation water and photoperiod both interact as limiting factors.
While more work is needed to disentangle the effects of the various climatic and environmental cues on greening trends, this study serves as a starting point for future research. For example, the dataset of disturbed pixels identified by the BFAST method, discarded here because beyond the scope of the paper, deserves a dedicated study addressing questions, such as the recovery speed after different types of disturbances, including fires, land-use changes, forest pest attacks, avalanches, landslides, etc. Additionally, as noted by [18], future studies on greening should expand the scale to an alpine-wide domain in order to: (i) confirm or discard patterns observed here and (ii) take into account processes operating at a larger scale, such as biogeographic gradients of precipitation, nitrogen deposition, and degree of disturbance. Finally, the combination of multiple sources of high-resolution satellite data can help addressing specific questions, such as impacts of drought and heat waves on greening trends or the dynamics of primary succession in the wake of glacier retreat.