Monitoring NDVI Inter-Annual Behavior in Mountain Areas of Mainland Spain (2001–2016)

: Currently, there exists growing evidence that warming is ampliﬁed with elevation resulting in rapid changes in temperature, humidity and water in mountainous areas. The latter might result in considerable damage to forest and agricultural land cover, affecting all the ecosystem services and the socio-economic development that these mountain areas provide. The Mediterranean mountains, moreover, which host a high diversity of natural species, are more vulnerable to global change than other European ecosystems. The protected areas of the mountain ranges of peninsular Spain could help preserve natural resources and landscapes, as well as promote scientiﬁc research and the sustainable development of local populations. The temporal statistical trends (2001–2016) of the MODIS13Q1 Normalized Difference Vegetation Index (NDVI) interannual dynamics are analyzed to explore whether the NDVI trends are found uniformly within the mountain ranges of mainland Spain (altitude > 1000 m), as well as in the protected or non-protected mountain areas. Second, to determine if there exists a statistical association between ﬁnding an NDVI trend and the speciﬁc mountain ranges, protected or unprotected areas are studied. Third, a possible association between cover types in pure pixels using CORINE (Co-ordination of Information on the Environment) land cover cartography is studied and land cover changes between 2000 and 2006 and between 2006 and 2012 are calculated for each mountainous area. Higher areas are observed to have more positive NDVI trends than negative in mountain areas located in mainland Spain during the 2001–2016 period. The growing of vegetation, therefore, was greater than its decrease in the study area. Moreover, differences in the size of the area between growth and depletion of vegetation patterns along the different mountains are found. Notably, more negatives than expected are found, and fewer positives are found than anticipated in the mountains, such as the Cordillera Cant á brica (C.Cant.) or Montes de Murcia y Alicante (M.M.A). Quite the reverse happened in Pirineos (Pir.) and Montes de C á diz y M á laga (M.C.M.), among others. The statistical association between the trends found and the land cover types is also observed. The differences observed can be explained since the mountain ranges in this study are deﬁned by climate, land cover, human usage and, to a small degree, by land cover changes, but further detailed research is needed to get in-depth detailed conclusions. Conversely, it is found that, in protected mountain areas, a lower NDVI pixels trend than expected (>20%) occurs, whereas it is less than anticipated in unprotected mountain areas. This could be caused by management and the land cover type.


•
Analyzing if significantly temporal NDVI short-term trends (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) occurred in the mainland Spanish mountains. • Studying if the possible patterns indicate a growing or decreasing trend for the different ranges due to their different climatic and vegetative conditions. • Studying if there exists an association between NDVI significant trends and location in protected or unprotected mountain areas. • Analyzing a possible association between NDVI trends and land cover types from mainland Spain.  Table 1 shows the percentage of protected area within every mountain area studied, and the mountain name abbreviation. Figure 1 shows the map of the present paper's study area. It must be noted that the lines in blue match the 1000 m elevation level and the black dotted lines indicate the total area of each mountain range. The mountain relief in mainland Spain involves discontinuities and contrasts in climate. The most important factors in the definition of mountain climate are altitude, latitude, orientation, massif, and continental. Generally, for every 100 m increment in elevation, the temperature drops 0.5 • C [47][48][49], which means greater rainfall possibilities (including in the form of snow), increased frost and shrinkage of the vegetative period.

Study Area
Mountain climate territories are those located more than 1000 m above sea level and, for that reason, only pixels with an altitude higher than 1000 m (see Figure 1) were the focus of this study. The mean annual temperature in those territories is lower than 10 • C, with mild summers (no month exceeds 22 • C) and cold winters (months with the average temperature below 0 • C). Additionally, there is abundant total annual precipitation (higher than 1000 mm) [47][48][49].
The study area can be divided into three major latitudinal regions: north peninsular mountains, continental mountains, and south peninsular mountains. Three floors or levels within these three major areas can be distinguished: subalpine forest floor (1000-1500 m), alpine forest floor (1500-2500 m) and snow forest floor (>2500 m).
North peninsular mountains receive a total annual precipitation greater than 1800 mm, where the subalpine floor forms a continuous band from M.L.G to the C.Cat. The typical vegetation here comprises conifers (Mediterranean zone) and deciduous forest (Atlantic zone). The average annual temperature is about 10 • C, and total annual precipitation exceeds 1100 mm. The maximum rainfall changes from winter (west) to autumn (east). The alpine forest floor (1500-2500 m) loses the continuity of the previous level and appears in isolated cores in the west and in the Pir., in which it reaches its maximum development. The temperatures are much more severe than those of the sub-alpine forest floor: annual averages are below 6 • C and, in these areas, the total precipitation usually is higher than 1200 mm. These general characteristics conceal pronounced anomalies, however, due to the orientation of the relief or to the different influences of the Atlantic and the Mediterranean. The shaded areas on the Mediterranean side keep the ground covered by snow for seven months, two more than in the sun, and the Atlantic influence doubles the rainfall (2500 mm) and better distributes it throughout the year, and raises average temperatures 2-4 • C. The vegetation of the alpine forest floor comprises conifers and deciduous forestland spaces and, in the highest altitude areas, meadows and shrubs.  Table 1.
Mountain climate territories are those located more than 1000 m above sea level and, for that reason, only pixels with an altitude higher than 1000 m (see Figure 1) were the focus of this study. The mean annual temperature in those territories is lower than 10 °C, with mild summers (no month exceeds 22 °C) and cold winters (months with the average temperature below 0 °C). Additionally, there is abundant total annual precipitation (higher than 1000 mm) [47][48][49].
The study area can be divided into three major latitudinal regions: north peninsular mountains, continental mountains, and south peninsular mountains. Three floors or levels within these three major areas can be distinguished: subalpine forest floor (1000-1500 m), alpine forest floor (1500-2500 m) and snow forest floor (> 2500 m).
North peninsular mountains receive a total annual precipitation greater than 1800 mm, where the subalpine floor forms a continuous band from M.L.G to the C.Cat. The typical vegetation here comprises conifers (Mediterranean zone) and deciduous forest (Atlantic zone). The average annual temperature is about 10 °C, and total annual precipitation exceeds 1100 mm. The maximum rainfall changes from winter (west) to autumn (east). The alpine forest floor (1500-2500 m) loses the continuity of the previous level and appears in isolated cores in the west and in the Pir., in which it  Table 1.
The snow floor, above 2500 m, is reduced at the center of the Pir. The average annual temperatures are below 0 • C, the seasons are diluted before a winter that lasts more than ten months, and summer is a warm and short parenthesis. Precipitation reaches 3000 mm in some areas, where mosses or lichens can be found [47][48][49][50][51].
Continental mountains show total annual precipitation between 1000 and 1500 mm. Average yearly temperatures range between 6 and 10 • C. The winters (0 • C/−2 • C) are long and cold and the summers are warm (although always below 22 • C). The orientation of the continental mountain relief is zonal, which entails an outstanding dissymmetry between the north and south slopes. The former are cooler and cloudier than the latter, yet receives a higher volume of rainfall because the southwest winds are more humid than those from the northwest. Due to their altitude, interior reliefs are framed within the sub-alpine forest floor. Finally, south peninsular mountains also show a subtropical mountain climate. The most outstanding characteristics which these zones present are lower winter temperatures than those of the previous climates mentioned. Summers are also warmer, although no month reaches 22 • C, which is widely exceeded in the surrounding lowlands. The average precipitation ranges between 800 and 1000 mm, but this amount oscillates significantly depending on the altitude and, especially, the proximity to the sea and exposure to the southwesterly winds. The maximum rainfall occurs in winter, followed by autumn and spring. The subtropical mountain suffers a deep drought in summer, like the rest of the territory of the peninsula to the southeast. Here, coniferous and deciduous forests and xerophytic plants are found [47][48][49][50][51].
Within these mountainous zones, we have also studied the protected zones, namely, the national parks and biogeographical regions in the continental zone and from mountains higher than 1000 m. Those protected areas represent a considerable percentage of the total mountainous area in mainland Spain and are a vast reservoir of natural richness.

CORINE (Co-Ordination of Information on the Environment) Land Cover Maps
The CORINE (Co-ordination of Information on the Environment) Land Cover (CLC) is a set of European environmental landscape maps [52]. These maps are a free download from the Land Copernicus project website. The vector version consists of a dataset at a scale of 1:100,000 and a minimum mapping unit of 25 ha [52].
The occupancy map of the CORINE Land Cover for 2006 (version 18.5) containing three levels of depth, which includes 44 types of land cover, is used in this study because it is from the middle of the period 2001-2016. The map of changes between 2000 and 2006 and between 2006 and 2012 is utilized to calculate percentages in the mountainous areas.

The NDVI from the Applied MODIS Dataset Used
The NDVI is a spectral index calculated from the reflectance in the red (R) and near-infrared (NIR) bands as follows [53]: NDVI = (NIR − R) / (NIR + R). The causal relationship between NDVI and the fraction of photosynthetically active radiation (FAPAR) or the net primary production of plants [54] has been demonstrated empirically for different systems [33][34][35]. The NDVI dataset of mountains within mainland Spain for the dates between 2001 and 2016 has been used.
The 16-day composite MODIS VI product from Terra (MOD13Q1), which shows a 250-m spatial resolution, is used. This dataset was downloaded from the Earth Data NASA website or the area of study (four tiles). The MODIS dataset was first processed by fusing the four tiles of every date between 2001 and 2016 (368 dates) with the help of an iterative process in TerrSet software [55]. Three hundred and sixty-eight images belonged to the 368 dates (one image every 16 days from 2001 to 2016, both inclusive).
The MOD13Q1 product provides a unique NDVI value (calculated with bands 1 and 2 of the MODIS sensor) every 16 days per pixel, allowing the composition of a gap-free temporal series. The MOD13Q1 provides a 16-day composite NDVI giving quality information from the layer of 250 m in terms of pixel reliability (PR), and acquisition dates from the 250 m 16 composite days of the year (CDOY) [56]. The PR layer informs about the pixels' reliability. The CDOY layer provides the date each vegetation index value was acquired [56]. The MODIS team uses an algorithm which chooses the best available pixel value from all the acquisitions from the 16-day period. The criteria utilized to select the best pixel value is low clouds, low view angle and the highest NDVI [56]. The MOD13Q1 product provided NDVI values with an approximate frequency of 16 days and an indicator of reliability to assure goodness in the next procedure before use in the annual time-series process.

Smoothing Normalized Difference Vegetation Index (NDVI) Pixel Curves in Spain with TIMESAT Software (Filtering Noise NDVI Profiles)
Using the TIMESAT software, Jönsson and Eklundh [57] used data from all the pixels in the study area from each year to determine the shape of the final year's curve.
TIMESAT provides processing methods based on least-squares fitting to the upper envelope of the vegetation indexes data (here the NDVI). An asymmetric Gaussian filtering method [57] was performed for calculating the NDVI from all pixels that made up the study area. The Gaussian method is an ordinary least-squares method, where data are fit to model functions. More detail about the functions that follow the TIMESAT software can be seen in Jönsson and Eklundh [57].
The PR information associated with every NDVI MODIS curve belonging to every pixel that comprised the study zone was considered. TIMESAT provides a weighting mechanism in which some values in the time series can be more influential than others, therefore, higher weights for high-quality retrievals and low weights for low-quality recoveries were assigned. Points assigned with lower weights in the smoothing procedure also have a minor influence on the curve fit [57].
Following getting the final smoothing curves of NDVI per pixel, the average of every pixel that comprises the study area (n = 368) per year (2001-2016) were calculated. Then, those pixels in which the NDVI value was less than zero were removed because they belonged to rocks and bare floor, not vegetation.

Temporal Trend Analyses
The average NDVI annual value per pixel from the smoothing NDVI pixel curves were calculated and described in Section 2.4 of the present paper, followed by the average NDVI annual pixel value which was evaluated in the serial correlation. The serial correlation defines how neighboring values in a series depend on each other. This serial correlation is usually used to detect trends in ecological time series [58].
A Theil-Sen estimator was used to examine the trends within each NDVI pixel smoothing curve (time-series) [59]. Theil-Sen's slope is considered a potential replacement of ordinary least squares for linear regression that computes the median of the pair wise slopes of all data cases that meet the ranking criteria. Thus, Theil-Sen's slope is not sensitive to outliers and is better than ordinary least squares for liner regression since it requires low a priori information of the measurement errors [60]. Next, the sign of the NDVI series trend was observed through the slope of the Theil-Sen regression performed. The positive-sign Theil-Sen slopes were those considered growing trends while the negative-sign Theil-Sen slopes were those considered browning (decreasing) trends of NDVI.
Subsequent to getting the Theil-Sen slopes from all the pixels of the study area, a Mann-Kendall test was used to detect the significance of the greening or browning trends [61,62]. The Mann-Kendall test is a robust non-parametric trend detection method. A significant trend was assumed when the z-statistic was less than 0.1 (90% confidence interval).
The Theil-Sen regression and the Mann-Kendall tests were performed using the Clark Labs TerrSet Earth Trends Modeler (ETM) [55]. ETM is an integrated suite of tools for the analysis of image time series data associated with Earth observation remote sensing images. Using these measures, one raster layer with the Mann-Kendall significance values and another raster layer with the Theil-Sen slopes from all the pixels of the study area were produced. Accompanying the last-mentioned dataset, with the help of a Geographic Information System (GIS) software, a layer in which the pixels with positive significant NDVI trends (growing NDVI trends (code = 1)), the pixels with negative NDVI trends (browning trends (code = −1)) and the pixels with non-significant trends (those with a Mann-Kendall z-statistic less than 0.1) (code = 0) was created. This layer, with the three classes (1, −1 and 0), was combined with the layer that contained the information of each mountain range (see Figure 1 and Table 1), as well as the layer with the protected or non-protected mountain zones and with the digital terrain model (DTM) from mainland Spain. Thus, pixels could be counted with significant positive (growing) and negative (browning) NDVI trends and without statistical trends per different Sustainability 2018, 10, 4363 7 of 24 mountains and per protected and non-protected areas. The chi-squared test of independence was used to determine if there was a significant relationship, first between the trend and the mountain categorical variables ( Table 2 second row for more detail) and second, between the trend and the protection rather than the mountain categorical variable (Table 2 third row for more detail). The mountain variable had 13 categories which matched the 13 different mountains from the study area; protection had two groups that were the protected and non-protected zones in the mountains; and, the trend variable had three categories, positive (growing NDVI), negative (browning NDVI) and no significant trend of NDVI. chi-squared tests between the trend and protected variables per mountain range were applied ( Table 2, fourth to final row for more detail). The frequency of each group for one variable was compared across the types of the second nominal variable for every test performed. The data can be displayed in a contingency table where each row represents a category for one variable, and each column represents a group for the other variable.
Graphs of NDVI trend distribution along different altitudes by mountains were created to assess the possible relationship between the NDVI trend and the altitude.
The possible relationship between the variables trend (with positive, negative and non-trend categories) and the type of land cover also was studied. The CORINE Land Cover Map of 2006 (level 3) filtered to have only pure pixels of the categories (82% pureness) was used. Pixels of the CORINE Land Cover Map were almost pure or belonged to one category, making a buffer of 125 m (half the pixel size of MODIS data used) in the vector CORINE Land Cover Map. This method was used to get NDVI pixel purity belonging to at least 82% of the land cover. This form of having pure pixels was used by Arrogante-Funes et al. [63,64].
A chi-squared test between the trend and cover type in which the cover type was built up by the most representative land cover categories was conducted. Those categories were non-irrigated arable land (Non-Irri.), agro-forestry areas (Agro-For), broad-leaved forest (Broad-leav.), coniferous forest (Con. For), natural grasslands (Nat. Grass) and transitional woodland-shrub (Tran. Wood-shr.). The chi-squared test was performed in a general way (with the whole study area) and then, a chi-squared test with the same variables used in the general approach per each mountainous area was conducted. The chi-squared analyses included only those mountain-land cover combinations in which there were at least 30 pure pixels. For the chi-squared test we used R 3.5 version free software.
Finally  Figure 2 indicates how the distribution of the percentage of negative and positive significant trends from 2001-2016 at different altitudes in mainland Spain are different than the total of pixels that make up mainland Spain. The positive trend curve seems to be different than the negative one. It also is observable how the curve with the percentage of pixels that do not show trends is similar to the curve that represents the total of pixels that make up mainland Spain. Generally, these results suggest how the altitude might be related to having pixels with NDVI significant trends. Focusing on the elevation higher than 1000 m, the present article study area could show how the positive trends curve is similar to the total percentage curve and the non-trend curve, whereas, the negative one has percentages higher than those other curves already mentioned. This fact is presented overall between 1300-1900 m of altitude.

Temporal Normalized Difference Vegetation Index (NDVI) Trend along Different Altitudes
were at least 30 pure pixels. For the chi-squared test we used R 3.5 version free software.
Finally, land cover changes by mountain with the CORINE Land Cover map of changes available from 2000-2006 and 2006-2012 were observed, and a Spearman rank correlation test between the areas with significant NDVI trends between 2001-2016 and the areas with land cover changes between 2000-2006 and 2006-2012, respectively, were performed. Figure 2 indicates how the distribution of the percentage of negative and positive significant trends from 2001-2016 at different altitudes in mainland Spain are different than the total of pixels that make up mainland Spain. The positive trend curve seems to be different than the negative one. It also is observable how the curve with the percentage of pixels that do not show trends is similar to the curve that represents the total of pixels that make up mainland Spain. Generally, these results suggest how the altitude might be related to having pixels with NDVI significant trends. Focusing on the elevation higher than 1000 m, the present article study area could show how the positive trends curve is similar to the total percentage curve and the non-trend curve, whereas, the negative one has percentages higher than those other curves already mentioned. This fact is presented overall between 1300-1900 m of altitude.

Figure 2
Relative percentage of pixels (Y-axis) with a significant NDVI trend along altitude in meters (X-axis) in the whole of mainland Spain. It must be noted that the percentage shown on the Y-axis refers to the total number of pixels belonging to all mainland Spain with a negative NDVI trend (red), the non-NDVI trend (yellow), positive NDVI trend (green) and the total pixels (blue). Figure 3 represents the percentage of pixels with negative or positive NDVI trends, non-trends and, finally, the total of pixels in each altitude higher than 1000 m, the mountainous zones within mainland Spain. Figure 3 shows how the curves from the negative and positive NDVI trends are different than the curve of distribution, which represents all the pixels in the mountainous systems. It shows at a glance how the 13 mountainous zones (Figure 3a-m) present differences in the trend distribution curves at different altitudes. It seems remarkable that, in general, the green line evolutions are more similar to the blue lines, than to the red lines. This implies that the greening (growth) in the vegetation patterns are near to the expected (pixels), whereas the red (decrease) in the vegetation trends are more different than expected. Relative percentage of pixels (Y-axis) with a significant NDVI trend along altitude in meters (X-axis) in the whole of mainland Spain. It must be noted that the percentage shown on the Y-axis refers to the total number of pixels belonging to all mainland Spain with a negative NDVI trend (red), the non-NDVI trend (yellow), positive NDVI trend (green) and the total pixels (blue). Figure 3 represents the percentage of pixels with negative or positive NDVI trends, non-trends and, finally, the total of pixels in each altitude higher than 1000 m, the mountainous zones within mainland Spain. Figure 3 shows how the curves from the negative and positive NDVI trends are different than the curve of distribution, which represents all the pixels in the mountainous systems. It shows at a glance how the 13 mountainous zones (Figure 3a-m) present differences in the trend distribution curves at different altitudes. It seems remarkable that, in general, the green line evolutions are more similar to the blue lines, than to the red lines. This implies that the greening (growth) in the vegetation patterns are near to the expected (pixels), whereas the red (decrease) in the vegetation trends are more different than expected.    Figure 4 shows how the positive NDVI trends of the pixels were more common than the negative patterns from 2001 to 2016. The percentage of pixels with significant negative trends were approximately 6.7%, whereas the portion of pixels with significant positive trends reached 11.4%. Sustainability 2018, 10, x FOR PEER REVIEW 10 of 24 Figure 4 shows how the positive NDVI trends of the pixels were more common than the negative patterns from 2001 to 2016. The percentage of pixels with significant negative trends were approximately 6.7%, whereas the portion of pixels with significant positive trends reached 11.4%. The I.S. with 26% and the S.C with 15.5% mountainous zones were those that had the most area from the total without significant NDVI trends in the study area. The C.Cant., with 29.37%, was the mountainous range with the most pixels from the total with a significant NDVI negative trend in the study area. Finally, the mountain systems with the greatest number of pixels with a significant positive trend were the I.N., with 19.8%, and the Pir. ranges, with 19.6%. These patterns can be seen in Figure 2b. Figure 2 demonstrates how there were mountains that had a predominance of significant NDVI negative trends, such as C.Cant. or L.G., whereas, it can be seen how I.N or M.V had a predominance of positive trends rather than negative ones. Another remarkable fact found in Figure 4 was that C.Cant., M.L.G. and I.S. mountains had most of the negative NDVI significant trends seem to match protected areas. Figures 5 and 6 are included to show the percentages of negative (red color bars) and positive (green color bars) trends, now calculated from the total area in each mountain ( Figure 5) and within every mountain from the total protected area (Figure 6a), and the percentage of negative and positive trend pixels from the total non-protected area. The key results can be summarized as follows: The I.S. with 26% and the S.C with 15.5% mountainous zones were those that had the most area from the total without significant NDVI trends in the study area. The C.Cant., with 29.37%, was the mountainous range with the most pixels from the total with a significant NDVI negative trend in the study area. Finally, the mountain systems with the greatest number of pixels with a significant positive trend were the I.N., with 19.8%, and the Pir. ranges, with 19.6%. These patterns can be seen in Figure 2b. Figure 2 demonstrates how there were mountains that had a predominance of significant NDVI negative trends, such as C.Cant. or L.G., whereas, it can be seen how I.N or M.V had a predominance of positive trends rather than negative ones. Another remarkable fact found in Figure 4 was that C.Cant., M.L.G. and I.S. mountains had most of the negative NDVI significant trends seem to match protected areas. Figures 5 and 6 are included to show the percentages of negative (red color bars) and positive (green color bars) trends, now calculated from the total area in each mountain ( Figure 5) and within every mountain from the total protected area (Figure 6a), and the percentage of negative and positive trend pixels from the total non-protected area. The key results can be summarized as follows: Sustainability 2018, 10, x FOR PEER REVIEW 11 of 24 Figure 5. Percentage of pixels (Y-axis) with significant NDVI positive (green) and negative (red) trends from the total pixels that make up every mountain range (X-axis). Contained in every positive and negative bar is the exact percentage of the NDVI trend from the total pixels in every protected and non-protected mountainous system. Abbreviations are reported in Table 1.  Table 1.

Temporal Normalized Difference Vegetation Index (NDVI) Trend along Different Mountains and Protected Areas
The mountainous system with the greatest percentage of significant NDVI trends from the total of each area was M.T.C.R. (more than 25% of all pixels of these mountains). By contrast, the mountainous systems with the lowest percentage of significant NDVI trend pixels were the I.S and the M.M.A.
It is also interesting to note that M.V, Pir., C.Cat., I.N and M.C.M. showed a great difference between the percentage of positive and negative pixels, with a predominance of positive ones (always Figure 5. Percentage of pixels (Y-axis) with significant NDVI positive (green) and negative (red) trends from the total pixels that make up every mountain range (X-axis). Contained in every positive and negative bar is the exact percentage of the NDVI trend from the total pixels in every protected and non-protected mountainous system. Abbreviations are reported in Table 1. Percentage of pixels (Y-axis) with significant NDVI positive (green) and negative (red) trends from the total pixels that make up every mountain range (X-axis). Contained in every positive and negative bar is the exact percentage of the NDVI trend from the total pixels in every protected and non-protected mountainous system. Abbreviations are reported in Table 1.  Table 1.
The mountainous system with the greatest percentage of significant NDVI trends from the total of each area was M.T.C.R. (more than 25% of all pixels of these mountains). By contrast, the mountainous systems with the lowest percentage of significant NDVI trend pixels were the I.S and the M.M.A.
It is also interesting to note that M.V, Pir., C.Cat., I.N and M.C.M. showed a great difference between the percentage of positive and negative pixels, with a predominance of positive ones (always Figure 6. Percentage of pixels (Y-axis) with significant NDVI positive (green) and negative (red) trends from the total pixels that make up every mountain range (X-axis). (a) The percentage is calculated from the total of pixels which built up the protected areas in every mountain. (b) The percentage is calculated from the total of pixels which built up the non-protected areas in every mountain. Contained in every positive and negative bar is the exact percentage of NDVI trends from the total pixels in each protected and non-protected mountainous system. Abbreviations are reported in Table 1.
The mountainous system with the greatest percentage of significant NDVI trends from the total of each area was M.T.C.R. (more than 25% of all pixels of these mountains). By contrast, the mountainous systems with the lowest percentage of significant NDVI trend pixels were the I.S and the M.M.A.
It is also interesting to note that M.V, Pir., C.Cat., I.N and M.C.M. showed a great difference between the percentage of positive and negative pixels, with a predominance of positive ones (always the percentage of positive trend pixels higher than 19). In contrast, the greatest difference between negative and positive pixels within a mountainous system occurred in C.Cant. and M.T.C.R. (see Figure 5), with more positive pixels.
It seems remarkable how in C.Cant., M.L.G and, to a lesser extent, in M.T.C.R., M.M.A. and B.S. mountain ranges, there was a higher percentage of negative NDVI trends (2001-2016) in protected areas than in non-protected. A lower rate of positive NDVI trends was found in protected areas compared to non-protected ones. Regarding M.V., C.Cat. and M.C.M., the reverse was the case-a higher percentage of positive trends in protected than in non-protected zones and a lower percentage of negative ones in protected compared with non-protected areas (see Figure 6) was found.
The results from the chi-squared tests (see Table 2) showed how the categorical variables Significant trend (positive, negative and non-trend) and Mountain range statistically were associated. It also associated the Significant trend with the Protection trend (protected or non-protected) when the test was performed considering all the pixels that comprise the study area and, also, when the test was made separately by mountainous systems. Having significant trends, and even the sign of the NDVI trend (growing or decreasing in the NDVI value), is related to the mountain studied and to being within a protected or non-protected area in these mountain areas. Figure 7 shows the graphs with the difference in percentage (X-axis) between the expected and the observed number of pixels with positive (growing), negative (decreasing) and no NDVI trend (equal). Figure 7a shows the differences by mountains and Figure 7b shows the differences calculated for the total mountainous area (>1000 m) including protected or not protected categories. Figure 8 shows the differences between expected and observed pixels with NDVI trends, but in every mountain-protected or mountain-not-protected zones.
Sustainability 2018, 10, x FOR PEER REVIEW 12 of 24 the percentage of positive trend pixels higher than 19). In contrast, the greatest difference between negative and positive pixels within a mountainous system occurred in C.Cant. and M.T.C.R. (see Figure 5), with more positive pixels. It seems remarkable how in C.Cant., M.L.G and, to a lesser extent, in M.T.C.R., M.M.A. and B.S. mountain ranges, there was a higher percentage of negative NDVI trends (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) in protected areas than in non-protected. A lower rate of positive NDVI trends was found in protected areas compared to non-protected ones. Regarding M.V., C.Cat. and M.C.M., the reverse was the case-a higher percentage of positive trends in protected than in non-protected zones and a lower percentage of negative ones in protected compared with non-protected areas (see Figure 6) was found.
The results from the chi-squared tests (see Table 2) showed how the categorical variables Significant trend (positive, negative and non-trend) and Mountain range statistically were associated. It also associated the Significant trend with the Protection trend (protected or non-protected) when the test was performed considering all the pixels that comprise the study area and, also, when the test was made separately by mountainous systems. Having significant trends, and even the sign of the NDVI trend (growing or decreasing in the NDVI value), is related to the mountain studied and to being within a protected or non-protected area in these mountain areas. Figure 7 shows the graphs with the difference in percentage (X-axis) between the expected and the observed number of pixels with positive (growing), negative (decreasing) and no NDVI trend (equal). Figure 7a shows the differences by mountains and Figure 7b shows the differences calculated for the total mountainous area (>1000 m) including protected or not protected categories. Figure 8 shows the differences between expected and observed pixels with NDVI trends, but in every mountain-protected or mountain-not-protected zones.   Table 1. Figure 7a shows how, from these mountainous systems that showed lesser observed pixels with NDVI positive trends in addition to higher observed pixels with NDVI negative trends, C.Cant., M.L.G. and M.T.C.R. were those in which the negative trends occurred more often than expected. Conversely, M.V., Pir, I.N., M.C.M., and C.Cat. showed fewer negative NDVI trends than expected, in addition to higher positive NDVI trends than expected by the statistics performed in the present paper (see also Figure 7a). Furthermore, in the cases of the last five mountain ranges mentioned, these differences between occurred and expected were higher than 45%. Figure 7b shows how, in protected areas, the percentage of the difference between observed and expected NDVI negative trend pixels in protected areas were higher than 20%, while, in non-protected  Figure 7a shows how, from these mountainous systems that showed lesser observed pixels with NDVI positive trends in addition to higher observed pixels with NDVI negative trends, C.Cant., M.L.G. and M.T.C.R. were those in which the negative trends occurred more often than expected. Conversely, M.V., Pir, I.N., M.C.M., and C.Cat. showed fewer negative NDVI trends than expected, in addition to higher positive NDVI trends than expected by the statistics performed in the present paper (see also Figure 7a). Furthermore, in the cases of the last five mountain ranges mentioned, these differences between occurred and expected were higher than 45%. Figure 7b shows how, in protected areas, the percentage of the difference between observed and expected NDVI negative trend pixels in protected areas were higher than 20%, while, in nonprotected regions the reverse was the case. Alternatively, the positive NDVI trends were very similar between the observed and expected counts. Figure 8a shows how, in the C.Cant. protected area, there were higher numbers of negative NDVI trend pixels than expected and fewer positive NDVI trend pixels than expected. The lastmentioned pattern occurred as a reverse case within the non-protected area of C.Cant. It happened in the same way for M.L.G; however, the percentage of the difference between expected and observed counts were lower than they were found within C.Cant. (See also Figure 8a,d). Figure 8i results reveal how the M.T.C.R in protected areas had a similar number of observed pixels to the NDVI trends to  Figure 8a shows how, in the C.Cant. protected area, there were higher numbers of negative NDVI trend pixels than expected and fewer positive NDVI trend pixels than expected. The last-mentioned pattern occurred as a reverse case within the non-protected area of C.Cant. It happened in the same way for M.L.G; however, the percentage of the difference between expected and observed counts were lower than they were found within C.Cant. (See also Figure 8a,d). Figure 8i results reveal how the M.T.C.R in protected areas had a similar number of observed pixels to the NDVI trends to those expected; however, in the non-protected areas, M.T.C.R showed a higher negative NDVI trend of observed pixels than expected and also lower positive ones than expected. Figure 8f,m showed how C.Cat. and M.C.M. had, in non-protected areas, more negative NDVI trend pixels than expected (difference between observed and expected counts higher than 50%) and, also, there was a lower percentage of positive than expected (difference between observed and expected counts higher than 40%). Regarding C.Cat., the difference between trend frequencies found and expected in the protected area were not high, nevertheless, in M.C.M. these differences were higher than in C.Cat. Considering the M.C.M. mountain range, the observed negative trends in the protected pixels were fewer than expected while, in contrast, the observed positive patterns were higher than expected.

Temporal Normalized Difference Vegetation Index (NDVI) Trends with Pure Pixels of Cover Types
The chi-squared results shown in Table 3 reveal how there also exists an association between the pure pixels with positive, negative or no trend and the land cover type. Table 4 shows the observed pixels by the most representative land covers in the study area. The most representative classes are the Con. For. and the Non-irri. followed by the Broad-leav. land cover type. The area of the total study area (3,084,612.94 ha in total) is observed with positive, negative and non-trend of NDVI by classes and the overall percentages of negative, positive and not trend are similar to those from the whole study area (as discussed in Section 3.2). The most common are the non-trend pixels, followed by the positive pixels which reached 10.16%, and then the negative trend pixels with a percentage of 7.34.  Figure 8 shows how Tran. Wood-shr. was the land cover type with the highest difference between observed and expected positive pixels (almost 80%). This difference revealed greater positive pixels than expected. Additionally, the number of negative pixels was lower than expected (difference smaller than 20%). Agro-for. and Non-Irri. were those with the highest difference between observed and expected NDVI negative trend pixels, with higher for anticipated than for observed. It also was found there were fewer positive than expected while Con. For., and the Broad-leav. were those cover types in which it was found the worst changes to the NDVI. Con. For. and the Broad-leav. showed greater negative NDVI trend pixels than expected. The NDVI evolution from 2001-2016 seems to be worse in Broad-leav. than in Con. For. because fewer positives were observed than expected while in Con. For greater positives were observed than expected (see Figure 9 for more detail).

Land Cover Change (2000-2006 and 2006-2012)
During the period between 2000 and 2006 changes in land cover in 0.89% of the total area were found, and from 2006 to 2012 changes in land cover were in 0.73% of the entire area (9,432,279.201 ha is the total area). Table 5 shows the percentage of change by mountains. It is observable how the mountains that had higher cover type changes were M.LG., I.S., S.C., and C.Cant. while the most changes were from forest categories to burnt areas and from regions burnt to shrub (these results are not shown in the text).  Table 1.

Discussion
This study analyzed short-term trends (growing and decreasing) of the NDVI over mountain areas (higher than 1000 m) in mainland Spain for the period from 2001 to 2016. The amount of space with growing (significant positive NDVI trends) and decreasing (significant negative NDVI trends) NDVI values were studied to see if differences were due to differences in the climatic, management and vegetation types of each mountain. Additionally, as the protection in mainland Spanish mountains is very important to preserve their vast biodiversity [16,[65][66][67], the study of the NDVI evolution, differentiating between protected and not protected areas, was included. Finally, the behavior of the NDVI during these years was explored for differences between cover types using pure pixels (one land cover category) and to see if the land cover change during the period studied could explain the patterns in the NDVI behavior.
First of all, in mainland Spain the NDVI trends seem to show different patterns at different altitudes. These differences also appeared between mountain ranges. The differences in the trends observed within different altitudes in the study zone are not a new find because, in the Mediterranean mountains and, generally, in other parts of the world, the vegetation succession processes and, thus, the vegetation distribution, are dependent on spatial and topographic patterns such as altitude or slope and orientation. Iterations between topographic effects on hydrological processes are some of the factors most widely used in differentiating vegetation patterns [68][69][70][71][72][73] for that reason.
Generally, in the studied mountain ranges (elevation higher than 1000 m) a predominance of positive significant NDVI trends rather than negative ones (11.4% versus 6.7%) was found. These percentages were similar when pure pixels of the most representative land cover types in the study area were used (see Table 5 for more detail). These results are in line with the global greening trends observed in other remote sensing vegetation studies around the world, [33,35,[74][75][76][77][78][79][80] and phenological field observations [81][82][83]. The studies cited provide evidence of an advance in the beginning of the growing season since three decades ago. The last study's trends are related to the increase in the temperatures of spring-summer and, in general, to the increase in the mean annual temperatures in mainland Spain and, in general, around the world. An accumulation of growing degree days concentrated, overall in the eastern United States and Europe, also were reported.
Regarding the percentages of positive and negative significant NDVI trends studied, differences were found between mountain ranges in the study area that are explained by the human, climatic and topographic (altitude) variables. These differences in the trends found also were related to the type of land cover in the mountain, as results showed. Generally, Tran. and Wood-shr. were found to be the classes with growing trends of NDVI while, by contrast, Coniferous Forest and Broad-leaved forest were those that suffered a decline in the NDVI patterns. Non-Irrigated arable land and Agro-forestry areas were found to be the most stable. This stability might be caused by the human influence that exists, compared to the other land covers. Differences for the same land cover between mountains indicate that the NDVI short-trends found in this article are explained, in part, by land cover types and, in part, by climatic factors (growth in mean temperatures, for instance [46]).
Additionally, changes in cover types from 2000 to 2006 were found in 0.89% of the total area, while, from 2006 to 2012, changes in land cover were found in 0.72% of the area. These changes were, overall, due to fires that had caused a change from the transitional shrub or forest categories to burnt areas or regions burnt to transitional shrub land covers; however, the magnitude is much lower than the magnitude of the NDVI trend found. Additionally, a statistical correlation between the area with change and the area with NDVI trends by mountains was not seen, so the cover type changes possibly might explain the trends. Nonetheless, the land cover change definitively was not the principal cause of those NDVI significant trends. Thus, the patterns are explained by cover type, land cover changes (to a lesser extent) and by climatic factors. These explanations also are supported by other authors who reported how vegetation dynamics are driven by climate, in particular by precipitation, incoming radiation and temperature [84][85][86][87][88][89], nutrient availability, and short-term natural and anthropogenic disturbances (fires, logging, insect epidemics, land-use change, and agricultural management, for example) that can be crucial at various spatiotemporal scales [90][91][92][93].
Specifically, in the Mediterranean mountains of Spain, the vegetation trends during the past few decades have been related to the abandonment of traditional activities [43,51,94,95]. Vidal-Macua et al. [95] discuss how land use management and the disturbances that are caused by the activities in the different land covers have been and continue to be the main driving force in vegetation succession. Conversely, in the Iberian Peninsula, there has been a remarkable increase in temperatures and a decrease in precipitation over the past few decades [96][97][98], which has led to an increase in the severity of droughts, especially in the Mediterranean area [99][100][101][102]. Higher quantities of decreasing vegetation were found in the Atlantic area than in the Mediterranean area, which was an unexpected pattern.
Looking at the statistical results, it is observable how there exists a significant relationship between finding negative, positive NDVI trends or not finding trends, from 2001 to 2016, and the mountain zone study. Generally, the positive occurred more frequently than expected based on the mean and the negative less frequently than was expected (see Figure 7a). These results match with other findings reported in the literature that discuss how positive vegetation activity is recorded mainly in mountain areas within Mainland Spain [43][44][45][46]103]. Papagiannopoulou et al. [104] showed how, for the period from 1981 to 2010, the production of vegetation in mainland Spain was controlled in the north by temperature and incoming irradiation and in the center, south, and east mostly by the precipitation regime. Thus, mountain zones have precipitation as their limiting factor. Karneli et al. [105] and Khorchani et al. [46] showed that, when water is the limiting factor for vegetation growth, the land surface temperature-NDVI correlation is expected to be negative, however, when the limiting factor is the energy, the relationship between the two variables is anticipated to be positive. Mountain areas in Spain are characterized mainly by high precipitation values and positive climatic water balances. Although these areas are affected by drought periods, overall, in spring and summer the vegetation continues developing because water availability is still enough; thus, slightly more significant positive trends than expected were found as well as less negative trends than expected. Some studies discuss how even arid or semi-arid areas have seen an increase in the photosynthetic activity due to the increase of CO 2 in the atmosphere as a result of human activities and the fertilization process [80,106].
It can be appreciated how, in mountain systems such as Montañas Vascas, Pirineos, Ibérico Norte, Cordillera Catalana, Sistema Central, and Béticas Sur and Montes de Cádiz y Málaga, a larger percentage of positive trends than expected was found and, in contrast, a smaller portion of negative trends than expected was found. These areas had seen increases in their photosynthetic activity during the time studied. These mountainous areas are within the Alpine and Mediterranean regions, so the increase in the temperature, joined with enough precipitation and tolerance to drought, resulted in more positive trends than expected. Mountains such as Cordillera Cantábrica, Montes de León y Galicia, Ibérico Sur, Montes de Toledo y Ciudad Real, Bética Sur o Montes de Murcia y Alicante had more pixels with negative trends and less with positive ones than expected. The photosynthetic activity began decreasing from 2001 to 2016 (see Figure 6a). The difference between the observable negative NDVI trends and the expected is notable, overall, in Cordillera Cantábrica, Montes de León y Galicia, Montes de Toledo y Ciudad Real y Montes de Murcia y Alicante.
Considering Cordillera Cantábrica, this pattern was observable mainly in the more Atlantic area, and the results are not expected because, in theory, the limiting factor was not the precipitation and the temperature has been slightly increasing in that area [46]. Currently, all research reported an increment in the NDVI value for the last three decades [42,79,107]. Khorchani et al. [46] showed how the correlation between NDVI and land surface temperature is not significant for most of the areas located in the northwest of Spain, so the temperature seems to not affect this pattern. The increase in the trees of the forest in the Cordillera Cantábrica (coniferous and broad-leaf forest) is controlled, in part, by the incoming energy and by the cloudy atmosphere. The climatic patterns in the northwest of Spain are controlled by the North Atlantic jet stream [108]. Recently, anomalies in the trajectory of the North Atlantic jet stream [109] were found that could be related to a decrease in the cloudy atmosphere or that could provoke more severe droughts in some critical months of the year for the growth of vegetation. Other authors have reported that Cordillera Cantábrica suffered a land degradation with high rates of land abandoned and fragmented, a severe overgrazing, and an increase in wildfire occurrence during the past few decades [110][111][112]. Other studies suggest that the Fagus sylvatica species in northern zones of Spain is less resistant to climatic change than those located in southern areas of Spain [113].
To conclude, all these factors explain finding more negative and less positive trends than expected in the northwest of Spain's mountain area. It was interesting to further the research in these areas, and not just in the semi-arid or arid zones for which the number of climatic change studies to date is higher than for humid areas. The results in M.L.G could be explained along with the results in Cordillera Cantábrica. However, land degradation is not as high and the species are in a range in which they need less humidity than those in Cordillera Cantábrica [113], a fact that could explain the percentage of negative NDVI trend pixels being lower than those found in the Cordillera Cantábrica. M.T.C.R and M.M.A are limited by precipitation and droughts. The increase in temperature and the drought suffered in the past few decades in these two mountain ranges could result in a decrease in vegetation activity. Additionally, in Montes de Murcia y Alicante, the abandonment of crop areas and urban growth has been blamed until now [79,107,114] and could explain, in part, these patterns.
To assess the relationship between the NDVI trends and the protected or non-protected areas, the study's statistical results show how there exists a significant association between the two variables. The negative NDVI trends occurred more frequently than expected in protected areas (with a difference between observed and expected higher than 20%) and, by contrast, less frequently than expected in non-protected areas (with a difference between observed and expected higher than 20%) (see Figure 7b for more information). These results match with those reported by Martinez-Fernández et al. [115], in which the authors observed a slight decline in natural surfaces in protected areas located in mainland Spain, with variations between the type of cover. They found that the transitional woodland shrub had a significant increase over 4% and, by contrast, more standing forest and open spaces with grasslands experienced a decline. This study's mountain protected areas have a higher percentage of standing timber and grasslands than the transitional woodland shrub and, maybe for that reason, the negative trends occurred with a higher frequency than expected and, on the contrary, occurred in the protected areas which had more standing forest and grasslands. Climatic and management factors could cause the finding of less vegetation activity than expected in the protected mountain areas during the past few decades [115], however, to have a greater level of understanding of the causes, it is necessary to conduct more specific research.
Concerning the protection effectiveness, Rodriguez-Rodriguez [17] and Martinez-Fernandez [116] reported recent results in which they showed high levels of effectiveness in protected Spanish sites (National and Natural Parks and Natura network 2000) since the development of urban and more natural cover in those protected areas. Both explained how more resources are needed for proper management and to prevent the effects of global change in peninsular Spain overall, and specifically for the Natura network 2000, considering the latest global financial crisis could affect the development and management of these protected areas during the period studied here. It was found in Catalonia specific natural management programs that could, in part, explain why more positive trends than negative were found overall in the Cordillera Catalana protected area (for more details see Figure 6).

Conclusions
This study showed the NDVI trends in the mountain area (elevation >1000 m) within mainland Spain, one of the world's most naturally biodiverse sites and also one of the most vulnerable to the effects of global change, using the MODIS 13Q1 dataset. There was a focus on protected and non-protected mountain areas because these areas can provide environmental sustainability in important mountain ranges. Additionally, the trends by cover type and the land cover changes that occurred from 2000 to 2006 and from 2006 to 2012 were analyzed to explain the trends found. The principal conclusions are summarized as follows.
Generally, a predominance of significant positive NDVI trends over negative ones was found, similar to the results reported for other European mountain systems, as well as Spain, during the last three decades. More positive and less negative trends than expected relative to the median were found. These results show how the trends observed were explained by land cover type and, to a small degree, by cover type changes (although the percentages of change are low). Different behaviors of the same land cover while located in different mountains indicated that climatic effects (temperature, and precipitation changes during the past few decades) and land cover changes reported in other studies caused the NDVI short-term trends in the current study area.
There were critical mountain ranges in which less positive or more negative than expected pixels were found. From past study results, Cordillera Cantábrica was expected to show a great majority of positive trends due to climatic conditions (increase of the land surface temperature and no higher episodes of annual precipitation accumulation). However, some aspects reported by other authors, such as anomalies in the North Atlantic jet stream, habitat fragmentations, land abandonment or the increase of the land surface temperature could lead to decreasing trends in vegetation that should be interesting for further study.
An association between NDVI trends in mountain protected areas or non-protected areas was found. Additionally, higher negative and less positive trends than expected were found in protected areas (with a difference higher than 20% between observed and expected). Less negative and higher positive than expected trend areas in non-protected sites (with a difference higher than 20% between observed and expected) were found. This general situation also occurred when each mountain was studied separately. The highest differences between observed and expected cases were recorded in Cordillera Cantábrica and Montes de Murcia y Alicante. These patterns could be due to the relationship with the land cover in these areas, and not mainly due to management issues. The present study characterizes the vegetation patterns occurring in the mainland Spanish mountains during the period 2001-2016 and their protected and not protected areas. It is necessary to further investigate, however, to understand and even predict what can happen in those essential reservoirs of service ecosystems in future years.