Assessing Vegetation Decline Due to Pollution from Solid Waste Management by a Multitemporal Remote Sensing Approach

Nowadays, the huge production of Municipal Solid Waste (MSW) is one of the most strongly felt environmental issues. Consequently, the European Union (EU) delivers laws and regulations for better waste management, identifying the essential requirements for waste disposal operations and the characteristics that make waste hazardous to human health and the environment. In Italy, environmental regulations define, among other things, the characteristics of sites to be classified as “potentially contaminated”. From this perspective, the Basilicata region is currently one of the Italian regions with the highest number of potentially polluted sites in proportion to the number of inhabitants. This research aimed to identify the possible effects of potentially toxic element (PTE) pollution due to waste disposal activities in three “potentially contaminated” sites in southern Italy. The area was affected by a release of inorganic pollutants with values over the thresholds ruled by national/European legislation. Potential physiological efficiency variations of vegetation were analyzed through the multitemporal processing of satellite images. Landsat 5 Thematic Mapper (TM) and Landsat 8 Operational Land Imager (OLI) images were used to calculate the trend in the Normalized Difference Vegetation Index (NDVI) over the years. The multitemporal trends were analyzed using the median of the non-parametric Theil–Sen estimator. Finally, the Mann–Kendall test was applied to evaluate trend significance featuring areas according to the contamination effects on investigated vegetation. The applied procedure led to the exclusion of significant effects on vegetation due to PTEs. Thus, waste disposal activities during previous years do not seem to have significantly affected vegetation around targeted sites.


Introduction
As early as the 1970s, a "new type of forest decline" [1] was observed almost simultaneously in various countries and in very different climatic and geo-pedological conditions and affecting different species of conifers and broad-leaved trees. The symptoms shown were not simply attributed to the classic biotic and abiotic stresses that cause a reduction in growth, a reduction in the leaf surface, the discoloration of the leaves, an arrest of or reduction in diameter increment, or a reduction in 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 vegetation functionality is now well known and established. Among the PTEs, heavy metals (HMs) play a fundamental role, so much so that the European Environment Agency (EEA) has highlighted the importance [2] of the continuous monitoring of HM emissions tions in plants [58]. Many studies have highlighted the ability of satellite remote sensing (RS) to identify geochemical stresses on plants induced by high HM concentrations through changes in the spectral responses in various ranges, in particular, visible and infrared (IR) [59][60][61]. In more detail, an increase in spectral reflectance in the red region and a significant reduction in the near infrared (NIR) was observed [62,63]. Thus, the use of Vegetation Indices (VI), based on the spectral transformation of two or more bands, tends to minimize the disturbances due to the reflectance contribution of soil, atmosphere and other external factors [64,65], highlighting the spectral aliquot 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 [69] and for the space-time comparison of photosynthetic activity [51,70]. In fact, NDVI is widely used for the evaluation of plant biomass, the leaf area index (LAI), photosynthetic activity and chlorophyll content; namely, all factors more or less strongly influenced by the presence of HMs.
Traditional methods for monitoring HM contamination involve soil and plant tissue analysis, providing very accurate data yet very onerous and expensive to sample; therefore, they are not applicable on a large scale [71,72]. On the other side, RS represents a valid alternative method for monitoring HM contamination [73][74][75]. It has been shown that HM contamination causes a reduction in the dry mass of plants even before alterations of the photosynthetic mechanisms occur [47]. Thus, RS techniques are potentially able to identify any contamination before clear evidence on vegetation [54]. Some RS studies are based on the implementation of models correlating the satellite spectral response with the variation of physiological and morphological parameters of vegetation as a result of different levels of contamination [76][77][78]. The main advantages of RS are to provide information in near real-time with high efficiency in identifying vegetation state changes and to be a non-destructive, low-cost methodology that provides large-scale data [79]. The negative consequences of high concentrations of HMs accumulated in plants can clearly manifest their effects over time [80]. Therefore, long-term observations of vegetation are of considerable importance [81]. In this perspective, a multitemporal analysis of RS data can be used to identify HM contamination in a wide area, further because the stresses induced by HM contamination have the characteristic of being stable in space and time [50]. In addition, the continuous, multitemporal and large-scale monitoring by RS allows policymakers to implement decision support systems for land management aimed at restoring ecosystem functionality and resilience capacity [4]. It should be noted that the effects induced by HM pollution could be confused with other factors (water stress, nutrient deficiency, diseases and infestations, climatic anomalies, mismanagement, etc.) affecting vegetation functionality and causing the same macroscopic effects [51,82,83] 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 years [50]. Starting from these assumptions, this work 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 variations due to external factors with respect to PTE 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 consisted of territorial areas affected by pollution issues, both in terms of overt contamination of environmental matrices as well as potential threats to the environment and public health. In particular, the following three sites (Figure 1) located in the north-western part of the Basilicata region (Italy) were analyzed:

•
The landfill of Aia dei Monaci, located in the municipality of Tito; • The landfill complex in the Montegrosso-Pallareta area of Potenza; • The former incinerator, later a waste transfer center, in Vallone Calabrese, Potenza.

Study Sites
The pilot sites for the multitemporal monitoring of vegetation consisted of territorial areas affected by pollution issues, both in terms of overt contamination of environmental matrices as well as potential threats to the environment and public health.
In particular, the following three sites ( Figure 1) located in the north-western part of the Basilicata region (Italy) were analyzed:

•
The landfill of Aia dei Monaci, located in the municipality of Tito; • The landfill complex in the Montegrosso-Pallareta area of Potenza; • The former incinerator, later a waste transfer center, in Vallone Calabrese, Potenza. The complex located in Aia dei Monaci was first used as a landfill for municipal solid waste (MSW) from 1994 to 2004, when the authorized volumes were exhausted, and then as a waste transfer station for MSW, serving the municipalities of the Bacino "Potenza Centro" from 2007 to 2014. At present, the transfer station is inactive, and no waste treatment or storage activities are carried out on the site. Its Site Characterization Plan is being implemented.
In the site of Montegrosso-Pallareta (Potenza), there is a complex of landfills whose work began in 1986 and that was responsible for waste disposal in May 1989. Since February 2009, it has housed an MSW transfer station inside of it, serving 18 municipalities in the hinterland of Potenza. Nowadays, the activities of the waste transfer station are suspended. The various basins are not in operation; therefore, they have been closed ac- The complex located in Aia dei Monaci was first used as a landfill for municipal solid waste (MSW) from 1994 to 2004, when the authorized volumes were exhausted, and then as a waste transfer station for MSW, serving the municipalities of the Bacino "Potenza Centro" from 2007 to 2014. At present, the transfer station is inactive, and no waste treatment or storage activities are carried out on the site. Its Site Characterization Plan is being implemented.
In the site of Montegrosso-Pallareta (Potenza), there is a complex of landfills whose work began in 1986 and that was responsible for waste disposal in May 1989. Since February 2009, it has housed an MSW transfer station inside of it, serving 18 municipalities in the hinterland of Potenza. Nowadays, 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 partial operation at the end of 2005 with the start of the testing procedures, the "hot tests" of the industrial plants, pending the completion of the authorization process for the exercise. The activity, up to that moment being 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 above-mentioned study areas are reported as "potentially polluted sites" in the "Registry of the sites subject to the remediation procedure" of the Basilicata region as analyses carried out over the years on various environmental matrices (groundwater, surface water, soil) have shown the following analytes to be out of range:  Table 1.

Satellite Data
The constellation of existing satellites and the characteristics of onboard sensors were analyzed to reconstruct the vegetation trends. The choice fell on the Landsat constellation as it is the longest earth observation (EO) program, providing data for the past 50 years, and the satellites of its missions have almost comparable sensors onboard, if applying appropriate intercalibration procedures. The good geometric resolution of Landsat images is able to identify with sufficient accuracy the changes occurring on the Earth's surface. Landsat 7 Enhanced Thematic Mapper Plus (ETM+) images were discarded from the analyses due to Scan Line Corrector (SLC)-Off causing data gaps since June 2003 that affect two of the three study areas almost entirely. Therefore, only Landsat 5 TM and Landsat 8 OLI images were considered (Table 2), acquired in the time span between 1990 and 2018 and referred to the summer months (generally July and August) when the percentage of cloud cover is near to zero 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 [84][85][86]. This algorithm, extensively tested to ensure its accuracy [87], has undergone several evolutions and in the current version includes the effects of polarization, which has a substantial role on method accuracy [88]. 6SV currently represents one of the most efficient algorithms for various multispectral sensors [89,90]. The subsequent cloud masking operation was necessary Remote Sens. 2022, 14, 428 7 of 22 in order to detect the presence of cloud cover and related shadows in order to identify the images to be excluded from the analysis in the study areas. The algorithm used was Fmask [91,92], which is particularly suitable for the correct identification of clouds and shadows in Landsat images [93,94]. Finally, the time series images were cropped on the whole study area, including all three sites of interest.

Analysis of the Vegetation Evolution
Starting from the Landsat data, a methodology for analyzing NDVI temporal trends was developed to obtain respective maps of the vegetation evolution/involution and maps of the environmental criticalities, both over the 1990-2018 timeframe ( Figure 3). The procedure adopted for the multitemporal analysis of processed Landsat data in the period 1990-2018 was the Theil-Sen median trend [95,96]. The non-parametric Theil-Sen estimator tends to produce accurate confidence intervals even when data are not normally distributed [97,98] and in the presence of heteroskedasticity [99,100]. The method is robust with regard to outliers [101,102], which are always present in a time series of RS data and 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 [103][104][105].
The procedure allows computing the slope of the line joining the medians of the trends calculated for each pair of points belonging to the dataset. Such a slope is calculated as follows [95,96]: where β is the slope between the data of two points in the time series, and xi and xj are 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 pixels over time. The final result provides the rate of change for the considered time step, e.g., year, month, week, etc. [106]. The Theil-Sen method tends to provide similar results to the Ordinary Least Squares (OLS) regression [107][108][109] for identifying the line slope when there are many observations (e.g., very long historical series and/or annual data on a daily, weekly basis, etc.) that exclude the effect of outliers [110][111][112]. For small observations, the median trend provides more robust results and is particularly suitable for the treatment of RS data, especially at medium and high spatial resolution and, therefore, with reduced observation frequency [113,114].
The main advantage of Theil-Sen is its "breakdown bound"; for a robust estimator, The procedure adopted for the multitemporal analysis of processed Landsat data in the period 1990-2018 was the Theil-Sen median trend [95,96]. The non-parametric Theil-Sen estimator tends to produce accurate confidence intervals even when data are not normally distributed [97,98] and in the presence of heteroskedasticity [99,100]. The method is robust with regard to outliers [101,102], which are always present in a time series of RS data and 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 [103][104][105].
The procedure allows computing the slope of the line joining the medians of the trends calculated for each pair of points belonging to the dataset. Such a slope is calculated as follows [95,96]: where β is the slope between the data of two points in the time series, and x i and x j are 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 pixels over time. The final result provides the rate of change for the considered time step, e.g., year, month, week, etc. [106]. The Theil-Sen method tends to provide similar results to the Ordinary Least Squares (OLS) regression [107][108][109] for identifying the line slope when there are many observations (e.g., very long historical series and/or annual data on a daily, weekly basis, etc.) that exclude the effect of outliers [110][111][112]. For small observations, the median trend provides more robust results and is particularly suitable for the treatment of RS data, especially at medium and high spatial resolution and, therefore, with reduced observation frequency [113,114].
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 by them is approximately 29% [106]. In such a case, the important implication is that the observations affected by strong inter-annual climatic oscillations or significant variations due to external factors are not taken into consideration compared to those to be observed.

Analysis of Environmental Criticalities
Once vegetation trends were identified, their significance was analyzed, identifying the areas that had 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 [115][116][117] is not only able to verify the existence of a trend but also to estimate its statistical significance (separating the hypothesis H0 = absence of trend from H1 = presence of monotonous trend). On the other hand, Z Mann-Kendall is not able to identify the trend monotonicity, unlike the Theil-Sen slope estimator. Instead, because the Theil-Sen estimator does not provide any information regarding trend significance, it is very often used in combination with the Mann-Kendall test [79,[118][119][120][121].
The Mann-Kendall Z test estimates the trend significance by measuring the magnitude of the relationship between two successive points. Assuming a time series t 1 , t 2 , . . . , t n , corresponding to a data series x 1 , x 2 , . . . , x n , the value S of the Mann-Kendall test [115,116] is: where n is the length of the time series, and x i and x j are the observations, respectively, at time i and j. When n ≥ 10, the Mann-Kendall test's statistical value S is similar to a normal distribution with a 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), where α is the level of significance, means that the time series shows a significant trend, increasing for S > 0 and decreasing for S < 0. The significance level considered for remotely sensed data is usually 0.05 [106,122,123]. For α = 0.05, Z(α/2) = 1.96. This means that for |Z| > 1.96, there is a significant trend of increase or decrease. The level of significance is particularly suitable for the analysis of a time series of RS data in order to verify the significance of NDVI trends [106,124,125].

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 each site's centroid. The consideration underlying this assumption was that any possible variation in the vegetation functional efficiency induced by the anthropogenic stresses was evident in an area not far from potentially polluting sources [126][127][128][129]. For a unit area of 900 m 2 (pixel size), the obtained maps identify the multitemporal variation of the vegetation efficiency at different levels of intensity and, thus, the environmental criticalities.

Maps of the Vegetation Evolution
The Theil-Sen slope map (Figure 4) was calculated for the NDVI, a vegetation index that is able to express any eco-physiological stresses on vegetation [42,130]. Therefore, starting from the NDVI Theil-Sen slope, a map of the vegetation evolution ( Figure 5) was created for each study area to identify any stresses and trace them to the factors that they induced. In more detail, this map on vegetation growth/degrowth rates in the period 1990-2018 identifies the following classes: • (pixel size), the obtained maps identify the multitemporal variation of the vegetation efficiency at different levels of intensity and, thus, the environmental criticalities.

Maps of the Vegetation Evolution
The Theil-Sen slope map ( Figure 4) was calculated for the NDVI, a vegetation index that is able to express any eco-physiological stresses on vegetation [42,130]. Therefore, starting from the NDVI Theil-Sen slope, a map of the vegetation evolution ( Figure 5) was created for each study area to identify any stresses and trace them to the factors that they induced. In more detail, this map on vegetation growth/degrowth rates in the period 1990-2018 identifies the following classes:    For the survey sites, consisting of circular areas with a radius of 1 km, the distribution percentages of the vegetation evolution classes are shown in Table 3.  For the survey sites, consisting of circular areas with a radius of 1 km, the distribution percentages of the vegetation evolution classes are shown in Table 3.
For Aia dei Monaci, the areas showing a decrease in vegetation were very small (about 2.5%), with consistent decreases affecting only six pixels (0.2%). Even the invariant areas were very limited in extent, while those that showed increases exceeded 96% in the timeframe considered, with very sustained increases reaching approximately 54%.
For the Montegrosso-Pallareta site, the total areas of decrease were approximately 15% and, therefore, they were certainly more consistent than in Aia dei Monaci, although the areas with a significant decrease were around 2%. Unlike how the Aia dei Monaci site falls within a distinct agro-forestry area characterized by clearly reduced land use changes, the Montegrosso-Pallareta site falls in a peri-urban area with a high complexity of landscape dynamics. In fact, since the 1990s, urbanization processes and land use changes have been faster and have involved larger surfaces. Thus, the agro-forestry land use dynamics were more sudden with an alternation of crop abandonment, cultivation of uncultivated areas, afforestation, etc. As expected, Vallone Calabrese acted similar to Montegrosso-Pallareta because the two sites are geographically very close and, therefore, are in similar microstational and territorial structure conditions. The vegetation decline affected just over 10% of the investigated land, while more than 86% of the surfaces showed a trend in positive evolution, of which over 53% were with very sustained growth rates.

Maps of Environmental Criticalities
The maps of environmental criticalities showing statistically significant decreases in the vegetation were elaborated through the Z Mann-Kendall methodology.
First of all, a regression analysis between the Z Mann-Kendall values and the Theil-Sen slope values [131] was performed in order to verify the validity of the procedure. The analysis was conducted on each site, considering a circular area with a radius of 1 km. Figure 6 shows a significant correlation between the variables under consideration for each site, with a determination index (R 2 ) always greater than 0.8, and even greater than 0.9 in the case of Montegrosso-Pallareta. In more detail, the functional model, R 2 , and the Standard Error of Estimate (SEE) of the linear regression between Z Mann-Kendall and Theil-Sen slope values for the three sites of interest are reported in Table 4. For Aia dei Monaci, the areas showing a decrease in vegetation were very small (about 2.5%), with consistent decreases affecting only six pixels (0.2%). Even the invariant areas were very limited in extent, while those that showed increases exceeded 96% in the timeframe considered, with very sustained increases reaching approximately 54%.
For the Montegrosso-Pallareta site, the total areas of decrease were approximately 15% and, therefore, they were certainly more consistent than in Aia dei Monaci, although the areas with a significant decrease were around 2%. Unlike how the Aia dei Monaci site falls within a distinct agro-forestry area characterized by clearly reduced land use changes, the Montegrosso-Pallareta site falls in a peri-urban area with a high complexity of landscape dynamics. In fact, since the 1990s, urbanization processes and land use changes have been faster and have involved larger surfaces. Thus, the agroforestry land use dynamics were more sudden with an alternation of crop abandonment, cultivation of uncultivated areas, afforestation, etc.
As expected, Vallone Calabrese acted similar to Montegrosso-Pallareta because the two sites are geographically very close and, therefore, are in similar microstational and territorial structure conditions. The vegetation decline affected just over 10% of the investigated land, while more than 86% of the surfaces showed a trend in positive evolution, of which over 53% were with very sustained growth rates.

Maps of Environmental Criticalities
The maps of environmental criticalities showing statistically significant decreases in the vegetation were elaborated through the Z Mann-Kendall methodology.
First of all, a regression analysis between the Z Mann-Kendall values and the Theil-Sen slope values [131] was performed in order to verify the validity of the procedure. The analysis was conducted on each site, considering a circular area with a radius of 1 km. Figure 6 shows a significant correlation between the variables under consideration for each site, with a determination index (R 2 ) always greater than 0.8, and even greater than 0.9 in the case of Montegrosso-Pallareta. In more detail, the functional model, R 2 , and the Standard Error of Estimate (SEE) of the linear regression between Z Mann-Kendall and Theil-Sen slope values for the three sites of interest are reported in Table 4.     are depicted for each study area in Figure 7. Further, some sampled points are also represented in these maps, and they are analyzed and discussed below.

Aia dei Monaci
Montegrosso-Pallareta Vallone Calabrese Functional Model Y = 0.000001 + 0.004313 * X Y = 0.000001 + 0.005762 * X Y = 0.000001 + 0.005951 * X R 2 0.82 0.87 0.94 SEE 0.00050 0.000318 0.000268 The detailed maps of environmental criticalities, evaluated through the Z Mann-Kendall test applied on the Theil-Sen median trend of NDVI in the time span considered (1990-2018) are depicted for each study area in Figure 7. Further, some sampled points are also represented in these maps, and they are analyzed and discussed below. Such maps identify the "non-sensitive areas", including portions of land that did not show any significant changes or showed a significant positive trend in the considered period, and "sensitive areas" in which, on the other hand, there was a statistically significant reduction in NDVI. For the three sites, the distribution of the areas of environmental criticality is shown in Table 5. For Aia dei Monaci, the negative variations, and therefore the areas of environmental criticality, were extremely low, being equal to 0.1% (only four pixels), while the areas without significant variations were approximately 23%. There were considerable areas (over 77%) with statistically significant positive variations. Analyzing the map (Figure  7a), it is possible to notice that the areas of environmental criticality were located in the northern part of the circular area of interest and were related to a quarry area already present at the beginning of the time series (1990) that has undergone a progressive expansion over the years.
In the case of the Montegrosso-Pallareta site, the largest surface (approximately 77%) fell into the "constant" class; unlike Aia dei Monaci, where the largest surfaces in Such maps identify the "non-sensitive areas", including portions of land that did not show any significant changes or showed a significant positive trend in the considered period, and "sensitive areas" in which, on the other hand, there was a statistically significant reduction in NDVI. For the three sites, the distribution of the areas of environmental criticality is shown in Table 5. For Aia dei Monaci, the negative variations, and therefore the areas of environmental criticality, were extremely low, being equal to 0.1% (only four pixels), while the areas without significant variations were approximately 23%. There were considerable areas (over 77%) with statistically significant positive variations. Analyzing the map (Figure 7a), it is possible to notice that the areas of environmental criticality were located in the northern part of the circular area of interest and were related to a quarry area already present at the beginning of the time series (1990) that has undergone a progressive expansion over the years.
In the case of the Montegrosso-Pallareta site, the largest surface (approximately 77%) fell into the "constant" class; unlike Aia dei Monaci, where the largest surfaces in the study area were affected by forest stands, in the site in question, the largest areas were 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 highlighted 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%) did not show a significant trend, while the critical areas (significant negative trend) were below 1%.

Discussion
In order to identify where the evolutionary trends of vegetation and, therefore, the environmental criticalities have occurred, both maps of the vegetation evolution and related environmental criticalities were intersected with the land use map of the Basilicata region ( Figure 2). The intersection of land use with the vegetation evolution classes for each study area is depicted in Figure 8. over the years.
For Vallone Calabrese, the percentage distribution of surfaces into criticality classes highlighted 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%) did not show a significant trend, while the critical areas (significant negative trend) were below 1%.

Discussion
In order to identify where the evolutionary trends of vegetation and, therefore, the environmental criticalities have occurred, both maps of the vegetation evolution and related environmental criticalities were intersected with the land use map of the Basilicata region ( Figure 2). The intersection of land use with the vegetation evolution classes for each study area is depicted in Figure 8.  Table 1.
For Aia dei Monaci (Figure 8a), especially moderate and strong vegetation decreases occurred very often in areas classified as "Non-irrigated arable land" (code 211). The cause was the decrease in NDVI due to a land use change from uncultivated lands (pasture and wooded pasture) to a resumption of cereal cultivation, which has notoriously lower NDVI average values. However, these areas represent only 0.7% of the entire area of Aia dei Monaci and, as will be seen below, show decreases that the statistical analysis considered insignificant. Further decreases were related to the expansion of urbanized areas, due in particular to the cultivation of a stone quarry located in the north of the study area. The most sustained increases, as expected, concerned areas with forest cover (stands of oak species). The land use class relating to "Transitional woodland-shrub" also showed markedly positive trends. These are areas in which the abandonment of crops or the reduction in livestock pressure has led to the expansion of forest areas. These are the typical areas of forest recolonization.
For Montegrosso-Pallareta (Figure 8b), the vegetation decreases basically affected two land use classes: "Non-irrigated arable land" (code 211) and "Coniferous forest" (code 312). In the first case, the decreases, which are not statistically significant, were due to the resumption of arable farming in previously uncultivated areas (pasture areas)  Table 1.
For Aia dei Monaci (Figure 8a), especially moderate and strong vegetation decreases occurred very often in areas classified as "Non-irrigated arable land" (code 211). The cause was the decrease in NDVI due to a land use change from uncultivated lands (pasture and wooded pasture) to a resumption of cereal cultivation, which has notoriously lower NDVI average values. However, these areas represent only 0.7% of the entire area of Aia dei Monaci and, as will be seen below, show decreases that the statistical analysis considered insignificant. Further decreases were related to the expansion of urbanized areas, due in particular to the cultivation of a stone quarry located in the north of the study area. The most sustained increases, as expected, concerned areas with forest cover (stands of oak species). The land use class relating to "Transitional woodland-shrub" also showed markedly positive trends. These are areas in which the abandonment of crops or the reduction in livestock pressure has led to the expansion of forest areas. These are the typical areas of forest recolonization.
For Montegrosso-Pallareta (Figure 8b), the vegetation decreases basically affected two land use classes: "Non-irrigated arable land" (code 211) and "Coniferous forest" (code 312). In the first case, the decreases, which are not statistically significant, were due to the resumption of arable farming in previously uncultivated areas (pasture areas) and the crop rotations between cereals and legumes. In the latter, the decrease was sometimes statistically significant.
Instead, forests (both broad-leaved trees and a part of conifer reforestation) and the "transitional woodland/shrub" were characterized by sustained increases. These are areas of forest recolonization, where the abandonment of cultivation has allowed the expansion of shrub and tree vegetation.
For Vallone Calabrese (Figure 8c), the decreases mainly affected arable land, for a similar reason as the Montegrosso-Pallareta site. The NDVI reduction in the "Broad-leaved forest" class (code 311), which was not statistically significant, is due to the management of some areas where this class is present. In particular, the riparian vegetation is subjected to thinning in order to maintain the efficiency of ditches and drains. The increases obviously affected the tree vegetation and also the recolonization stage and the "Natural grasslands" (code 321), i.e., the grazing areas where the reduction in anthropogenic pressure (crop abandonment, reduction in grazing pressure) has led to the increase in herbaceous biomass over the years. The intersection between land use classes and maps of environmental criticalities (Figure 9) allows identifying clearly whether vegetation perturbations have occurred over the years due to potentially polluting anthropic activities.
as of forest recolonization, where the abandonment of cultivation has allowed the expansion of shrub and tree vegetation.
For Vallone Calabrese (Figure 8c), the decreases mainly affected arable land, for a similar reason as the Montegrosso-Pallareta site. The NDVI reduction in the "Broadleaved forest" class (code 311), which was not statistically significant, is due to the management of some areas where this class is present. In particular, the riparian vegetation is subjected to thinning in order to maintain the efficiency of ditches and drains. The increases obviously affected the tree vegetation and also the recolonization stage and the "Natural grasslands" (code 321), i.e., the grazing areas where the reduction in anthropogenic pressure (crop abandonment, reduction in grazing pressure) has led to the increase in herbaceous biomass over the years. The intersection between land use classes and maps of environmental criticalities (Figure 9) allows identifying clearly whether vegetation perturbations have occurred over the years due to potentially polluting anthropic activities.  Table 1.
For Aia dei Monaci (Figure 9a), a statistically significant decrease occur in "Artificial surfaces", that is, as already discussed above, the areas located in the north of the site that are characterized by the expansion of a quarry cultivation. The positive changes (significantly increasing trend) are instead related to forest stands, in particular to broadleaved forests (in detail, oak trees) and to areas of forest recolonization ("Transitional woodland-shrub"). The invariant areas fundamentally belong to the "Non-irrigated arable land" class, where there is an annual cycle of cultivation on the same surfaces.
For the Montegrosso-Pallareta site (Figure 9b), the positive trends occur mainly in correspondence with natural surfaces: from forest stands (both broad-leaved and part of the reforested areas with conifers) to the areas of forest recolonization following the reduction in anthropogenic pressure (crop abandonment) and natural herbaceous vegetation areas (pastures). The invariant areas involve both the surfaces with crops (nonirrigated arable land) and the small areas without any vegetation cover. Some conifers' reforestations show an arrest in growth, which had already begun to decline since 1995. A more careful examination of these reforestations allows identifying that the constant decline in terms of biomass is not related to the possible effect due to waste disposal ac-  Table 1.
For Aia dei Monaci (Figure 9a), a statistically significant decrease occur in "Artificial surfaces", that is, as already discussed above, the areas located in the north of the site that are characterized by the expansion of a quarry cultivation. The positive changes (significantly increasing trend) are instead related to forest stands, in particular to broadleaved forests (in detail, oak trees) and to areas of forest recolonization ("Transitional woodland-shrub"). The invariant areas fundamentally belong to the "Non-irrigated arable land" class, where there is an annual cycle of cultivation on the same surfaces.
For the Montegrosso-Pallareta site (Figure 9b), the positive trends occur mainly in correspondence with natural surfaces: from forest stands (both broad-leaved and part of the reforested areas with conifers) to the areas of forest recolonization following the reduction in anthropogenic pressure (crop abandonment) and natural herbaceous vegetation areas (pastures). The invariant areas involve both the surfaces with crops (non-irrigated arable land) and the small areas without any vegetation cover. Some conifers' reforestations show an arrest in growth, which had already begun to decline since 1995. A more careful examination of these reforestations allows identifying that the constant decline in terms of biomass is not related to the possible effect due to waste disposal activity, but to the partial failure of these stands due to local conditions, in particular pedological characteristics, and to the inadequate choice of plant species. This is, in fact, a common situation in the reforestation of Basilicata, where the choice of alien species that are not very suitable for extreme climate and above all soil features has led to the constant deterioration of reforested areas over the years, up to the partial or total failure of the plantation.
In the Vallone Calabrese site (Figure 9c), the invariant areas, similarly to the Montegrosso-Pallareta site, consist of arable land and pasture areas, while the positive trends basically occur both in forest areas and forest recolonization areas. Instead, the critical areas, statistically significant negative trends, are the "Artificial surfaces".
The analysis of NDVI Theil-Sen multitemporal trends of some points of interest (Figures 6 and 10) clarifies the vegetation dynamics in each investigated area.  In Aia dei Monaci, the points P9 and P10 fall into critical areas and show decreasing trends. The initial NDVI values in 1990 are already very low (about 0.1) and typical of rocky areas with sparse herbaceous vegetation. These values tend to constantly decrease due to the expansion of quarry cultivation over the years. No other areas of environmental criticality can be identified on this site. Then, there are areas (P5) characterized by a constant trend over the years and points that, although showing a slightly increasing (P4 and P7) or decreasing trend (P8), do not appear to have a statistical significance with regard to the Z Mann-Kendall test. As mentioned above, most areas have a statistically significant positive variation in the NDVI (P2 and P6) and are attributable to forest stands.
In Montegrosso-Pallareta, the positive trends are related to forests (P12) also in the recolonization stage (P14) or areas, such as the landfill body itself (P11), characterized by the presence of herbaceous and shrub artificial species, as a result of restoration operations (capping and greening). Points between P18 and P20 correspond to the reforestation stands affected by a forest fire in 2015, as evidenced by a significant decrease in NDVI. However, the growth decline started in 1995 because these plants, as already clarified above, are located in very difficult local conditions, characterized by heavily clayey asphyxiated soils on very high slopes.
For Vallone Calabrese, the invariant points (P22 and P24) consist of arable land and pasture areas, while the positive trends are mainly related to forests, forest recolonization areas (P26 and P25) and pastures (P21) where the reduction in anthropogenic pressure (grazing) has led to a biomass increase. The critical areas (significantly negative trends) are attributable to urbanization and infrastructural works (P28) and, partially, to arable land as a consequence of the restarting of the cultivation of pasture or wooded pastures on small surfaces with consequent biomass reduction (P29).
In summary, the analysis carried out on the three sites in order to identify the presence of the effects on the vegetation of any pollution leads to the conclusion that a significant alteration of plant systems' functionality with a consequent decline in the vegetation were not found in the areas of interest during the period under examination. The significant vegetation decline observable in limited areas of each site was due to landscapes dynamics and the consequent land use changes, or external disturbances with respect to the phenomenon under investigation.
The proposed methodology, which involves the application of the Theil-Sen slope and Z Mann-Kendall test to verify the statistical significance of trends, has proven to be particularly suitable in the processing of collected satellite data and taking into account the perturbations induced on vegetation by different drivers than those to be analyzed. In fact, over a long period of time and with few observations, NDVI variations may be due to external causes (disturbances) able to significantly affect its trends. For instance, a forest fire occurring near Montegrosso-Pallareta in 2015 (Figure 11a) or forest exploitation at Aia dei Monaci (Figure 11b), due to the usual forest management, represent outliers with respect to vegetation trends. In such cases, the Theil-Sen median trend is not affected, unlike OLS regression, by these extreme events, which represent aberrant values compared to the targeted phenomenon.
In fact, over a long period of time and with few observations, NDVI variations may be due to external causes (disturbances) able to significantly affect its trends. For instance, a forest fire occurring near Montegrosso-Pallareta in 2015 (Figure 11a) or forest exploitation at Aia dei Monaci (Figure 11b), due to the usual forest management, represent outliers with respect to vegetation trends. In such cases, the Theil-Sen median trend is not affected, unlike OLS regression, by these extreme events, which represent aberrant values compared to the targeted phenomenon.

Conclusions
The aim of this work was to verify the possibility of analyzing the effects induced on plant systems in areas surrounding sites where potentially polluting activities take place or have been carried out. The analyses concerned sites located in southern Italy (Basilicata region) intended for the disposal of municipal solid waste and where PTE concentrations above the threshold values permitted by current national legislation were recorded at some stage of their management.
The proposed analyses identified the multitemporal NDVI trends and their significance factually and efficiently, thus allowing the elaboration of maps of the vegetation evolution and related environmental criticalities. The methodology made it possible to evaluate the vegetation evolution per unit area from 1990 to 2018. Ultimately, the adopted workflow was able to discover the presence or absence over time of stressors related to anthropogenic waste disposal activities that could have influenced the significant alteration of plant systems' functionality. The applied methodology was able to exclude negative effects on the vegetation due to potentially polluting activities. In fact, the environmental criticalities were not very spatially large near the three sites of interest. In more detail, in Aia dei Monaci, the critical areas were almost absent, representing about 0.1% (0.36 ha) of the entire study area (ab. 314 ha), while in Montegrosso-Pallareta and Vallone Calabrese, the areas with a significant decrease in vegetation were 0.7% (2.0 ha)

Conclusions
The aim of this work was to verify the possibility of analyzing the effects induced on plant systems in areas surrounding sites where potentially polluting activities take place or have been carried out. The analyses concerned sites located in southern Italy (Basilicata region) intended for the disposal of municipal solid waste and where PTE concentrations above the threshold values permitted by current national legislation were recorded at some stage of their management.
The proposed analyses identified the multitemporal NDVI trends and their significance factually and efficiently, thus allowing the elaboration of maps of the vegetation evolution and related environmental criticalities. The methodology made it possible to evaluate the vegetation evolution per unit area from 1990 to 2018. Ultimately, the adopted workflow was able to discover the presence or absence over time of stressors related to anthropogenic waste disposal activities that could have influenced the significant alteration of plant systems' functionality. The applied methodology was able to exclude negative effects on the vegetation due to potentially polluting activities. In fact, the environmental criticalities were not very spatially large near the three sites of interest. In more detail, in Aia dei Monaci, the critical areas were almost absent, representing about 0.1% (0.36 ha) of the entire study area (ab. 314 ha), while in Montegrosso-Pallareta and Vallone Calabrese, the areas with a significant decrease in vegetation were 0.7% (2.0 ha) and 0.3% (1.0 ha), respectively. As widely discussed, such critical areas can be traced back to external factors (disturbances) rather than to pollution due to waste disposal activities.
In conclusion, it should be emphasized that the methodological approach adopted in this work has a general validity; it is not site-specific, thus it is exportable in other territorial contexts. It can be suitably applied to other activities with a high impact on the environment (e.g., mining, oil extraction, exploitation of groundwater, etc.) or natural phenomena that may induce eco-physiological stress on vegetation, affecting its functional efficiency and dynamics. Namely, the impacts that can be induced on ecological systems by global change in relation to phenology or productivity, such as forest fires and subsequent vegetation recovery, the risk of desertification, urbanization processes and all kinds of man-made pollution that can affect ecosystem functionalities.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are openly available here: http: //www.webgis-simbiosi.it/, accessed on 13 November 2021.