Assessing vegetation decline due to pollution from waste cycle activities by a multitemporal remote-sensing approach

: The work consisted in identifying possible effects from heavy metals (HMs) pollution due to waste disposal activities in three potentially polluted sites located in Basilicata (Italy), where a release of pollutants with values over the thresholds allowed by the Italian legislation was detected. The potential variations in the physiological efficiency of vegetation have been analyzed through the multitemporal processing of satellite images. In detail, Landsat 5 Thematic Mapper (TM) and Landsat 8 Operational Land Imager (OLI) images were used to calculate the Normalized Difference Vegetation Index (NDVI) trend over the years. Then, the multitemporal trends were analyzed using the median of Theil-Sen, a non-parametric estimator particularly suitable for the treatment of remote sensing data, being able to minimize the outlier effects due to exogenous factors. Finally, the subsequent application of the Mann-Kendall test on the trends identified by Theil-Sen slope allowed the evaluation of trends significance and, therefore, the areas characterized by the effects of contamination on vegetation. The application of the procedure to the three survey sites led to the exclusion of the presence of significant effects of HMs contamination on the vegetation surrounding the sites during the years of waste disposal activities.


Introduction
As early as in the 1970s, a "new type of forest decline" [1] was observed in various countries almost simultaneously, in very different climatic and geo-pedological conditions and affecting different species of conifers and broad-leaved trees. The symptoms shown were not simply attributable to the classic stresses of a biotic and abiotic nature and related to reduction of growth, reduction of the leaf surface, discoloration of the leaves, arrest or reduction of diameter increment, reduction of root phytomass and crown transparency. The phenomenon, whose causes have not yet been fully established, is still the subject of scientific research but the contribution of Potential Toxic Elements (PTEs) in the decline of the functionality of vegetation is now well known and established. Among the PTEs, heavy metals (HMs) play a fundamental role, so much as the European Environment Agency (EEA) has highlighted the importance [2] of a continuous monitoring of HM emissions in the environment at European level. Such monitoring is a priority of international programs, like the International Cooperative Program on Assessment and Monitoring of Air Pollution Effects on Forests (ICP Forests) and the International Cooperative Program on Integrated Monitoring of Air Pollution Effects on Ecosystems (ICP IM) of the United Nations Economic Commission for Europe (UNECE) and there are now numerous similar programs present at national level in various parts of the world [3,4].
Although HM pollution began with the industrial revolution [5], it is in recent decades that there has been an increase in the concentration of HMs in the soil that is causing serious environmental damage as a result of the problems of contamination. The quality of the environment and the activities of living organisms, in fact, are influenced by contamination from HMs [6][7][8]. Where the excessive concentration of these metals concerns environments accessible to animals or plant species of food interest, animal and human health can also be compromised [9][10][11].
HMs can be released into the environment both for natural drivers and as a result of anthropogenic activities. The main causes of emissions are anthropogenic sources, in particular mining activities [12][13][14], industry and energy production [15,16], the agricultural sector (fertilizers, pesticides, etc.), the use of fossil fuels [17], waste disposal [18][19][20], population growth and concentration in large urban areas [21,22]. The demographic and industrial increase has led, in fact, to the huge production of waste that is treated, stored, collected and disposed, involving risks for the environment in the various phases [23,24]. The controlled or uncontrolled disposal of waste can cause, due to accidental causes, the leakage and migration of contaminants, including inorganic ones, such as heavy metals, with a significant impact on the environment [25]. In some cases, the metals emitted continue to persist in the environment, for many years, even after the anthropogenic activities are decommissioned.
HMs are absorbed by plants, together with essential elements, from the soil, giving rise to potentially risky bioaccumulation phenomena not only for flora and fauna but also for human health. The toxicity of HMs in plants depends on the plant species, the chemical element, its concentration and the chemical-physical characteristics of the soil. The elements that can be significantly accumulated in plants are Cr, Mn, Fe, Ni, Cu, Zn, Cd, Pb and Hg. The problems caused by HMs are well known on plants [26][27][28][29][30]: growth reduction [31][32][33][34] until reaching high levels of toxicity capable of causing death [35]. High concentrations of HMs, in fact, damage plant structures, also through the replacement of essential ions with free radicals that affect the protein and enzymatic structure [36]. Furthermore, high concentrations of HMs in the soil negatively interfere with the radical absorption of essential nutrients, affecting physiological processes, such as respiration and photosynthesis. HMs also favor the production of Reactive Oxygen Species (ROS) by self-oxidation [37] which, in turn, results in negative effects on lipids and proteins of the cell membrane, on permeability, on Deoxyribo-Nucleic Acid (DNA), on photosynthetic properties ultimately causing physiological stress, alterations in phenology and growth reduction [38][39][40]. Further negative effects are attributable to physiological processes such as germination, accumulation and remobilization of reserve substances during germination [41].
The functional alterations of plants subjected to HMs pollution result in morphological alterations concerning the content of the pigments [42], the structure of the mesophyll, the characteristics of the leaf surface [43,44] and canopy architecture [45]. From the physiological point of view, the variations concern the transpiration rate and photosynthesis [46,47], stomatal conductance [48,49] and consequently the leaf temperature [50,51]. This essentially translates into effects ranging from the reduction in the number of leaves, to the reduction of the leaf surface and green pigments, and therefore to chlorosis [52], necrosis, leaf epynasty, red-brown discoloration [3,4,53]. Chlorosis is probably the most common phenomenon in response to heavy metal pollution stresses [33], causing a reduction in chlorophyll absorption around the spectral ranges of 680 nm and 550 nm [54].
Such a plant response affects the reflectance spectra of foliage and canopy under an increasing geochemical stress [55,56]. However, the changes in terms of reflectance can be species-dependent [57] and be influenced by environmental variables that affect the concentration of HMs in plants [58]. Numerous studies have highlighted, however, the ability of satellite remote sensing (RS) to identify geochemical stresses on plants induced by high concentrations of HMs, through changes in the spectral responses in various ranges, visible and infrared (IR) in particular, of the electromagnetic spectrum [59][60][61]. The studies showed, in particular, an increase in spectral reflectance in the red region and a significant reduction in the near infrared (NIR) [62,63]. The use of Vegetation Indices (VI), based on the spectral transformation of two or more bands, tends to minimize the reflectance disturbances due to the reflectance contribution of soil, atmosphere and other external factors [64,65], highlighting the spectral contribution of vegetation [66][67][68].
The NDVI Vegetation Index, based on the normalized difference of the Red and NIR bands, which are more sensitive to the variations induced by contamination, is particularly suitable for the multitemporal identification of the physiological efficiency of vegetation and for the space-time comparison of photosynthetic activity [51,69]. NDVI is in fact widely used for the evaluation of plant biomass, leaf area index (LAI), photosynthetic activity and chlorophyll content, all factors that are more or less strongly influenced by the presence of HMs.
Traditional methods for monitoring HMs contamination, therefore, involve soil and plant tissue analysis, and provide very precise data but are nevertheless onerous and expensive and not applicable on a large scale [70,71]. RS represents a valid alternative method for monitoring HMs contamination [72][73][74]. In addition, it has been shown that the HMs contamination causes a reduction in the dry mass of plants even before alterations of the photosynthetic mechanisms occur [47] and therefore the ability of RS to identify any contamination before that is visually evident on vegetation [54]. Studies using RS are based on the construction of models that correlate the spectral response of the satellite with the variation of physiological and morphological parameters of vegetation as a result of different levels of contamination [75][76][77]. Ultimately, the advantage of RS is to provide information in near real-time, with high efficiency in identifying changes in the state of the vegetation, and to be a non-destructive, low-cost methodology that provides largescale data [78]. The negative consequences of high concentrations of HMs, which result in their accumulation in plants, can clearly manifest their effects over time [79]. Therefore, long-term observations of vegetation are of considerable importance [80]. In this perspective, multitemporal analysis of RS data can be used to identify HMs contamination in a given large area, also due to the fact that the stresses induced by the contamination of HMs have the characteristic of being stable in space and time [50]. In addition, the continuous, multitemporal and large-scale monitoring allowed by RS allows policy-makers to implement more quickly tools for the management and governance of the territory aimed at restoring the functionality of the ecosystem and its resilience capacity [4]. Finally, it should be noted that the effects induced by HM pollution could be confused with other factors (water stress, nutrient deficiency stress, diseases and infestations, climatic anomalies, mismanagement, etc.) that affect the functionality of the vegetation and which can have the same macroscopic effects [51,81] and therefore very similar spectral variations. However, the latter stressors have a limited effect over time while the toxicity due to HM bioaccumulation persists over the years [50]. This work started from these assumptions and led to: a) multitemporal analysis of vegetation change in areas surrounding potentially polluted sites, through the study of NDVI trends, b) identification of a statistical procedure for analyzing the physiological trends of vegetation that did not take into account the variations due to external factors with respect to HMs pollution, c) analysis of the statistical significance of the multitemporal trends of NDVI for the possible identification of areas of environmental criticality due to the effect of contamination.

Study sites
The pilot sites for the multitemporal monitoring of vegetation consist of territorial areas affected by high environmental criticalities, both in terms of overt contamination of environmental matrices, as well as potential threats to the environment and public health. In the site of Montegrosso-Pallareta (Potenza) there is a complex of landfills whose works began in 1986 and which was responsible for the waste disposal in May 1989. Since February 2009 it has housed a MSW transfer station inside of it serving 18 municipalities in the hinterland of Potenza. As of today, the activities of the waste transfer station are suspended. The various basins are not in operation; therefore, they have been closed according to the procedures established by the regulations in force at the time of filling the basins.
The former incinerator in Vallone Calabrese (Potenza), built between 1988 and 2003, only came into operation partially at the end of 2005 with the start of the testing procedures with the "hot tests" of the industrial plants, pending the completion of the authorization process for the exercise. The activity, up to that moment provisional and never fully operational, ended in 2007. Today, the structures that house the industrial plants are not operational and there is no waste management and / or handling. The area is used as a depot for snow vehicles and other equipment for the ordinary operation of the "Azienda Comunale per la Tutela Ambientale" (ACTA), the local Municipal Environmental Protection Agency.
The three sites are reported as "sites potentially polluted" in the database of the Basilicata Region [82], as analyses conducted in various times and on various matrices (groundwater, surface water, soil) have shown the presence of analytes out of range, especially HMs (iron, nickel, lead, aluminum, selenium, etc.).

Satellite data
To reconstruct the vegetation trends, the constellation of existing satellites and the characteristics of the sensors on board were analyzed. The choice fell on the Landsat constellation as it was the only one that has been providing earth observation (EO) data for the past 50 years, with the various satellites that have almost comparable sensors on board, with the appropriate intercalibrations. The good geometric resolution of the images provided, among other things, is able to identify with sufficient accuracy the changes occurring on the Earth's surface. The careful analysis of the dataset of the various Landsat missions considered led to the exclusion of Landsat 7 Enhanced Thematic Mapper Plus (ETM+) images, as due to Scan Line Corrector (SLC)-Off since June 2003 there are data gaps that affect two of the three study areas almost entirely. The choice therefore fell on Landsat 5 TM and Landsat 8 OLI datasets (Table 1), starting from data from 1990 and identifying images in the summer period (generally July and August) that did not present, to the visual analysis, cloud cover on survey sites. The atmospheric correction of the images, absolutely necessary when dealing with satellite image time series, was carried out using the 6SV algorithm [83][84][85]. The algorithm, extensively tested to ensure its accuracy [86], has undergone several evolutions and in the current version includes the effects of polarization which has a substantial role on the accuracy of the method [87]. 6SV currently represents one of the most efficient algorithms for various multispectral sensors [88,89]. The subsequent cloud masking operation was necessary in order to verify whether, in the areas of interest, there was the presence of cloud cover and related shadows that would have led to the exclusion of the image from the analysis. The algorithm used is Fmask [90,91], particularly suitable for the correct identification of clouds and shadows of Landsat images [92,93]. Finally, the image stack of the time series was cropped on the study area in such a way as to include all the three sites of interest.

Analysis of the evolutionary trends of vegetation
Starting from the Landsat data, a methodology for analyzing the temporal trends of the NDVI vegetation index was developed in order to obtain respectively maps of the evolution / involution of vegetation and maps of the environmental criticalities, both over the 1990-2018 timeframe.
The procedure used for the multitemporal analysis of processed Landsat data in the period 1990-2018 is the trend of the Theil-Sen median [94,95]. The non-parametric Theil-Sen estimator tends to produce accurate confidence intervals even when data are not normally distributed [96] and in presence of heteroskedasticity [97,98].
Furthermore, the method is robust with regard to outliers [99,100] which are always present in time-series of RS data and which represent a "disturbance" with respect to the analysis of vegetation trends. Furthermore, in RS this non-parametric procedure is considered robust with respect to the seasonality, non-normality and frequently present autocorrelation both at intra and interannual scale [101][102][103].
The procedure calculates the slope of the line joining the medians of the trends calculated for each pair of points belonging to the dataset. From an analytical point of view, the slope is calculated as follows: where β is the slope between the data of two points in the time-series and xi and xj the corresponding values between the two points i and j (for j < i). The procedure is applied independently to each pixel and calculates the median of the slopes between all the n (n-1) / 2 pair combinations of the pixels over time. The final result provides the rate of change for the considered time step, e.g., year, month, week, etc. [104]. The Theil-Sen method tends to provide similar results to the Ordinary Least Squares (OLS) procedure [105][106][107], for identifying the slope of the line, when there are many observations (very long historical series and annual data of a daily, weekly type, etc.) as the effect of the outlier [108,109]. For small observations, the trend of the median provides more robust results and is particularly suitable for the treatment of RS data, especially at medium and high resolution and, therefore, with reduced temporal resolution [110,111].
The main advantage of Theil-Sen is its "breakdown bound": for a robust estimator the number of outliers that can be present in the dataset before the trend is affected is approximately 29% [104]. In such a case, the important implication is that the observations affected by strong inter-annual climatic variations or significant variations due to external factors are not taken into consideration compared to those to be observed.

Analysis of environmental criticalities
Once the trends (increases or decreases) of the vegetation have been identified, their significance was analyzed, i.e., the areas that have undergone environmental stress over the years such as to induce a significant reduction in physiological efficiency.
To carry out this verification, the Mann-Kendall test was used starting from the Theil-Sen slope trends. In fact, the Z Mann-Kendall test [112][113][114] is not only able to verify the existence of a trend, but also to identify its significance. Since the Theil-Sen estimator does not provide any information regarding the significance of the estimated trend, the nonparametric Mann-Kendall test is often used [115,116] associated with it. On the other hand, Z Mann Kendall despite being able to estimate the significance of a trend (separating the hypothesis H0 = absence of trend from H1 = presence of monotonous trend) is not able to identify by himself, if an increasing or decreasing trend is manifesting: thus the Theil-Sen slope (able to identify increases / decreases) is very often used with the Mann-Kendall test (able to identify the significance of trend) [78,117,118].
The Mann-Kendall Z test is able to measure the significance of the trend by measuring the magnitude of the relationship between two successive points. Assuming a time series t1, t2, …, tn, corresponding to a series of data x1, x2, …, xn, the value S of the Mann-Kendall test is: where n is the length of the time series, and xi and xj are the observations respectively at time i and j. When n ≥ 10, the Mann-Kendall test statistical value S is similar to a normal distribution with mean equal to 0. The variance of S is: The Mann-Kendall Z value is used to identify where the trend is significant: |Z| > Z(α/2) means that the time series shows a significant trend, increasing for S greater than 0 and decreasing for S less than 0; α is the level of significance. For a significance level of 0.05, Z(α/2) = 1.96. This means that for |Z| > 1.96 there is a significant trend of increase or decrease. The significance level considered for remotely sensed data is usually 0.05 [104,119,120]. The level of significance is particularly suitable for the analysis of timeseries of RS data in order to verify the significance of NDVI trends [104,121].

Results
The analyses for the identification of the multitemporal trends of the vegetation were conducted on a circular area with a radius of 1 km starting from the site's centroid. The consideration underlying this assumption was that any variations in the functional efficiency of vegetation induced by the stresses of the anthropic activities under study were, if present, evident in an area not far from potentially polluting sites [122][123][124]. Ultimately, the maps obtained identify, for a unit area of 900 m 2 (pixel size), the multitemporal variation of the efficiency of vegetation functionality at different levels of intensity and, finally, of the environmental criticalities.

Map of the evolutionary trends of vegetation
The Theil-Sen slope map (Figure 2) was developed for NDVI, a vegetation index able to express any eco-physiological stresses of vegetation [42,125]. Therefore, starting from the NDVI Theil-Sen slope, a map of the evolution of the vegetation which expresses the growth / decrease rates of the vegetation in the period 1990-2018 was obtained for each study area in order to identify any stresses and trace them to the factors that they induced. This map has been elaborated on the entire study area comprising the survey sites ( Figure  3) and identifies respectively the involution of the vegetation in 3 classes (slight, moderate and strong decrease) and the evolution of the same in 3 classes (slight, moderate and strong increase). Finally, the intermediate class, i.e., that of the invariant areas, is classified as "constant".  For the survey sites, consisting in circular areas with a radius of 1 km, the distribution percentages of the vegetation evolution classes are shown in Table 2. For Aia dei Monaci, the areas showing a decrease in vegetation are very small (about 2.5%) with consistent decreases affecting only 6 pixels (0.2%). Even the invariant areas are very limited in extent, while those that show increases exceed 96% in the time frame considered, with very sustained increases reaching approx. 54%.
For the Montegrosso-Pallareta site, the total areas of decrease are approx. 15% and, therefore, they are certainly more consistent than in Aia dei Monaci, although the areas with a significant decrease are around 2%. In fact, while the Aia dei Monaci site falls within a distinctly agro-forestry area, with decidedly contained dynamics of land use, in the case of the Pallareta site we are faced with a peri-urban area with a complexity of strongly elevated landscape with sustained evolutionary dynamics. In this case, in fact, the works attributable to urbanization and "land use" have been, since the 1990s, faster and have involved larger surfaces. The dynamics of use of the agro-forestry land were, consequently, more sudden with an alternation of crop abandonment, cultivation of uncultivated areas, afforestation, etc.
As we can legitimately expect, the Vallone Calabrese situation is similar to the Montegrosso-Pallareta site, because the two sites are geographically very close and, therefore, in similar microstational and territorial structure conditions. The declining areas affect just over 10% of the investigated land while more than 85% of the surfaces show a trend in positive evolution, with over 53% of the surface with very sustained growth rates.

Environmental criticality map
The maps of the environmental criticalities showing the statistically significant decreases of the vegetation were made with the Z Mann-Kendall methodology.
To verify the validity of the procedure, a regression analysis was performed between the Z Mann-Kendall values and the Theil-Sen slope ones [126]. The analysis was conducted on individual sites, namely circular areas with a radius of 1 km. Figure 4 shows for each site a significant correlation between the variables under consideration, with a determination index (R 2 ) always greater than 0.8 and even greater than 0.9 in the case of the Pallareta site. Functional model, R 2 and Standard Error of Estimate (SEE) of the regression between Z Mann Kendall vs. Theil-Sen slope values for the three sites are reported in Table 3.  The environmental criticality maps, using the Z Mann-Kendall test on the median trend series of NDVI Theil-Sen slope in the time span considered , are depicted in Figure 5, where some points of particular importance are also highlighted, which will be discussed below. The maps identify the "non-sensitive areas" (divided into portions of land that did not show any significant changes in the period considered and those with a significant positive trend) and "sensitive areas" in which, on the other hand, there is a statistically significant reduction in NDVI. For the three sites, the distribution of the areas of environmental criticality is shown in Table 4. For Aia dei Monaci, the negative variations, and therefore the areas of environmental criticality, are extremely low, being equal to 0.1% (only 4 pixels), while the areas without significant variations are approx. 23%. There are considerable areas (over 77%) with statistically significant positive variations. Analyzing the map (Figure 5a), it is possible to notice that the areas of environmental criticality are located in the northern part of the circular area of interest and are related to a quarry area already present at the beginning of the time-series (1990) and that has undergone a progressive expansion over the years.
In the case of the Montegrosso-Pallareta site, the largest surface (approx. 77%) falls into the "constant" class: unlike Aia dei Monaci where the largest surfaces in the study area are affected by forest stands, in the site in question, the largest areas are occupied by herbaceous crops (in particular arable land) which are renewed annually and, therefore, do not show significant variations in terms of quantity of biomass over the years.
For Vallone Calabrese, the percentage distribution of surfaces into criticality classes highlights a situation that is sufficiently similar to the Montegrosso-Pallareta site, as the structure of the landscape is similar in terms of anthropogenic, cultural and natural components. The larger areas (over 60%) do not show a significant trend while the critical areas (significant negative trend) are below 1%.

Discussion
In order to identify where the evolutionary trends of vegetation and, therefore, environmental criticalities have occurred, the maps of evolution / involution and environmental criticalities with land use were cross compared, using the Land Use map of Basilicata Region. The classes of land use that fall within the sites of interest are shown in Table 5. The intersection of land use with the classes of evolution / involution of vegetation is illustrated in Figure 6. For Aia dei Monaci (Figure 6a), the decreases (especially moderate and strong) occurred in areas currently in "Non-irrigated arable land" (code 211). The cause is attributable to a decrease in NDVI due to a change in land use: that is, there has been a shift from areas of neglected cultivation, attributable to the typology of "pasture" and "tree pasture", to a resumption of cereal cultivation, which has notoriously average values lower than NDVI. However, these areas represent only 0.7% of the entire area of Aia dei Monaci and, as it will be seen below, show decreases that the statistical analysis considers insignificant. Further decreases relate to the expansion of urbanized areas, attributable in particular to the cultivation of a stone quarry located in the north of the study area. The most sustained increases, as expected, concern areas with forest cover (stands of oak species). The land use class relating to "Evolving arboreal and shrub vegetation" also shows markedly positive trends: these are areas in which the abandonment of crops or the reduction of livestock pressure has led to the expansion of forest areas. These are the typical areas of forest recolonization.
For Pallareta (Figure 6b) we emphasize that the decreases basically affect two classes of land use: "Non-irrigated arable land" (code 211) and "Coniferous forest" (code 312). In the first case, the decreases, which are not statistically significant, are due to the resumption of arable farming in previously uncultivated areas (pasture areas) and the alternation of crop rotations between grasses and legumes. For the "Coniferous woods", however, in some cases the decrease is statistically significant. The topic will be explored later when we talk about the Maps of environmental criticalities.
With regard to the areas in which sustained increases are recorded, in addition to woods (both broad-leaved trees and part of conifer reforestation), the areas with "transitional woodland/shrub" should be noted: these are areas of forest recolonization, where the abandonment of cultivation has allowed the expansion of shrub and arboreal vegetation.
For Vallone Calabrese (Figure 6c) we remark that the decreases mainly affect arable land, for a similar argument of the Montegrosso-Pallareta site. The reduction in the "Broad-leaved forest" class (code 311), which is not statistically significant, is attributable to the management of some areas where this crop typology is present, in particular to the management of riparian vegetation, subjected to cleaning operations and thinning to maintain the efficiency of ditches and drains. The increases obviously affect the arboreal vegetation, also in the recolonization phase, and "Natural grasslands" (code 321), i.e., the grazing areas where the reduction of anthropogenic pressure (crop retreat, reduction of grazing) has led to the increase in years of herbaceous biomass. The intersection between land use data and the maps of environmental criticalities (Figure 7) makes it possible to identify with greater clarity whether damage to vegetation has occurred over the years due to the exercise of potentially polluting anthropic activities. For Aia dei Monaci (Figure 7a), the areas that have undergone a statistically significant decrease are those attributable to "Artificial surfaces", that is, as already discussed above, to the areas located in the north of the site, characterized by the expansion of the cultivation of a quarry. The positive changes (significantly increasing trend) are instead attributable to forest stands, in particular to broad-leaved forests (oak wood, in the specific case) and to areas of forest recolonization ("Transitional woodland-shrub"). The invariant areas are fundamentally attributable to "Non-irrigated arable land", where there is an annual cycle of cultivation on the same surfaces.
For the Montegrosso-Pallareta site (Figure 7b), the positive trends are mainly attributable to the "natural" surfaces: forest stands (both broad-leaved and part of the reforested areas with conifers), to the areas of forest recolonization following the reduction of anthropogenic pressure (crop abandonment) and areas with natural herbaceous vegetation (pastures). The invariant areas involve both the surfaces with herbaceous crops (non-irri-gated arable land) and the small areas without any vegetation cover. Even some reforestations of conifers show a stagnation in the growth, which, in some cases, as we will see later, becomes even a decline as early as the year 1995. A more careful examination of these reforestations has allowed us to identify that the constant decline, in terms of biomass, both to be related not to the possible effect due to the waste disposal activity, but to the partial failure of these plants due to the microstational conditions, pedological conditions in particular, and to the inadequate choice of plant species. This is, in fact, a common situation in reforestation in Basilicata where the choice of alien species, not very suitable for extreme climatic and, above all, pedological conditions, has led to constant deterioration over the years up to the partial or total failure of the plant.
In the Vallone Calabrese site (Figure 7c), the invariant areas, similarly to the Pallareta site, are attributable to arable land and pasture areas while the positive trends are mainly attributable to forest areas and forest recolonization areas. The critical areas, statistically significant negative trends, are attributable to the "Artificial surfaces'' class.
The analysis of NDVI's Theil-Sen multitemporal trends of some points of particular interest (Figures 5 and 8) clarifies the dynamics of the evolution / involution of vegetation.
In Aia dei Monaci, the points that fall into critical areas (points P9 and P10) are those that show decreasing trends. The initial NDVI values (dating back to 1990) are already very low (NDVI equal to about 0.1), a typical value of rocky areas with sparse herbaceous vegetation, and which over the years tends to constantly decrease due to the expansion of quarry cultivation. In the area of the Aia dei Monaci site, no other areas of environmental criticality can be identified. There are areas (point P5) in which there is no trend (neither positive nor negative) over the years and points that, although showing a slight upward trend (points P4 and P7) or downward (point P8), do not appear to have a statistically significant trend with regard to the Z Mann-Kendall test. As mentioned above, most areas have a statistically significant positive variation in the NDVI index (points P2 and P6) and are attributable to forest stands.
In Montegrosso-Pallareta, the positive trends are related to the woods (point P2) also in the recolonization phase (point P23) or to areas, such as the body of the landfill itself (point P21), which were affected, following the restoration operations (capping and greening), from artificial herbaceous and shrub species. Points P28 and P32 correspond to the reforestations which were affected by a forest fire in 2005, as evidenced with a significant decrease in NDVI. However, the growth decline started in 1995, because these plants, as already clarified above, are located in very difficult microstational conditions, characterized by heavily clayey asphyxiated soils on very high slopes.
For Vallone Calabrese, the invariant points (points P22 and P24) are attributable to arable land and pasture areas while the positive trends are mainly attributable to forest areas, forest recolonization areas (P16 and P15) and pastures (P11) in which the reduction of anthropogenic pressure (grazing) has led to a biomass increase. The critical areas (significant negative trends) are attributable to urbanization and infrastructural works (point P18) and, partially, to arable land as a consequence of the cultivation of grazing areas or tree pastures on small surfaces with consequent biomass reduction (point P19). In summary, the analysis carried out on the three sites, in order to identify the presence of effects on the vegetation of any pollution, leads to the conclusion that, in the areas of interest, physiological stresses with consequent decline of the vegetation are not found during the period under examination. The significant negative variations found in some points, in the various sites, are attributable to the dynamics of the landscape and to the consequent variations in land use or to external disturbances with respect to the phenomenon investigated.
The procedure used, which involves the application of Theil-Sen slope and Z Mann Kendall test to verify the statistical significance of the trends, has proved to be particularly suitable in the treatment of data collected by satellite and in taking into account, over a period of time relatively long, the variations induced on the vegetation by different drivers than those to be investigated. In fact, over a long period of time and with sufficiently reduced observations, the variations of NDVI may be due to causes external to the phenomenon investigated and such as to significantly affect it. In fact, as it can be observed (Figure 9), in the case of a fire occurred in 2005 at the Montegrosso-Pallareta site ( Figure  9a) or forest utilization cuts, as in the Aia dei Monaci site (Figure 9b) due to normal forest management, and which therefore represent outliers with respect to the analysis of vegetation trends, the Theil-Sen median trend procedure is not affected while the OLS tends to be strongly influenced by these extreme events which represent aberrant values compared to the phenomenon investigated.

Conclusions
The analyses conducted allow the identification of multitemporal trends and the measurement of their significance in a very punctual and efficient manner, allowing to obtain maps of the evolution of vegetation and environmental criticalities. The methodology made it possible to identify, pixel by pixel, what happened, from 1990 to 2018, on these unitary surfaces in relation to the evolution of the vegetation. Ultimately, the procedure made it possible to identify the presence or absence over time of stressors related to anthropogenic waste disposal activities that in some way could have influenced the physiological efficiency of the vegetation and, therefore, its functional stresses. The methodology, applied on the survey sites, allowed to exclude negative effects on the vegetation due to potential concentrations of pollutants. The areas of environmental criticality are greatly reduced in the survey sites and attributable, as widely discussed, to external factors (disturbances) with respect to possible pollution due to waste disposal activities.
In conclusion, it should be emphasized that the procedure used in this work is a general procedure, which is exportable in other territorial areas and for other activities that are presumed to have an influence on the vegetation by inducing stress, or for the evaluation of large-scale phenomena which can affect the functional efficiency of vegetation and its dynamics. Namely, the variations that can be induced on ecological systems by Global Change in relation to phenological variations [127,128] or in productivity [81,129,130], the phenomenon of forest fires and of subsequent recovery processes [131][132][133][134][135][136], the risk of desertification [137], urbanization processes [138,139] up to all those pollution phenomena of various types [140] that can induce variations in functionality of the ecosystems. Funding: The present work was a further result of the project "SIMBioSi" (2017-2019), co-financed by the Basilicata Region. The APC was funded by Centro di Geomorfologia Integrata per l'Area del Mediterraneo.

Data Availability Statement:
The data presented in this study are openly available here: http://www.webgis-simbiosi.it/.