Robust Satellite-Based Identiﬁcation and Monitoring of Forests Having Undergone Climate-Change-Related Stress

: Climate-induced drought events are responsible for forest decline and mortality in different areas of the world. Forest response to drought stress periods may be different, in time and space, depending on vegetation type and local factors. Stress analysis may be carried out by using ﬁeld methods, but the use of remote sensing may be needed to highlight the effects of climate-change-induced phenomena at a larger spatial and temporal scale. In this context, satellite-based analyses are presented in this work to evaluate the drought effects during the 2000s and the possible climatological forcing over oak forests in Southern Italy. To this aim, two approaches based on the well-known Normalized Difference Vegetation Index (NDVI) were used: one based on NDVI values, averaged over selected decaying and non-decaying forests; another based on the Robust Satellite Techniques (RST). The analysis of the ﬁrst approach mainly gave us overall information about 1984–2011 rising NDVI trends, despite a general decrease around the 2000s. The second, more reﬁned approach was able to highlight a different drought stress impact over decaying and non-decaying forests. The combined use of the RST-based approach, Landsat satellite data, and Google Earth Engine (GEE) platform allowed us to identify in space domain and monitor over time signiﬁcant oak forest changes and climate-driven effects (e.g., in 2001) from the local to the Basilicata region scale. By this way, the decaying status of the Gorgoglione forest was highlighted two years before the ﬁrst visual ﬁeld evidence (e.g., dryness of apical branches, bark detachment, root rot disease). The RST exportability to different satellite sensors and vegetation types, the availability of suitable satellite data, and the potential of GEE suggest the possibility of long-term monitoring of forest health, from the local to the global scale, to provide useful information to different end-user classes.


Introduction
Recent studies (e.g., [1][2][3]) report an increase in tree mortality due to climate change. Climate-induced tree mortality is affecting not only semi-arid regions but also mesic forests. It is a global phenomenon that impacts monsoonal savannas as well as conifer and broadleaves forests with a Mediterranean and tropical climate [4]. More frequent, long-lasting, severe events of drought and heat stress are among climate-driven causes of forest decline and die-off. Drought and/or heat-induced forest mortality is documented in several cases, in the last decades, all around the world [4].
Deadly effects on vegetation may be evident even after years from drought stress periods [4], but stress-induced tree mortality is usually preceded by functional and structural changes (see, for instance, growth changes [5]). Vegetation response to drought can be Land 2022, 11, 825 2 of 18 different on a local scale because several factors, such as soils, topographic position, tree density, and root systems, may play a different role [6,7]. Drought and heat stress may lead to a leaf area reduction and crown dieback up to a broad-scale forest die-off [8]. In the western Mediterranean region, recent studies [9] indicate drought as the main driver of the oak decline. Apart from site-specific conditions, some oak species are more drought-tolerant than others. Some authors [10] suggest, for example, the following increasing tolerance rank: Quercus robur (less tolerant), Quercus cerris, Quercus frainetto, Quercus pubescens, Quercus ilex, Quercus suber (more tolerant). A decline is documented in different Italian regions both in the north (Piedmont [11]) and in the south (in Sicily [12] and Basilicata [9,13]). In particular, in Southern Italy, several oak species have shown decline and mortality since the 2000s due to drought stress [10].
Stress analysis may be carried out by using more or less complex field methods such as the computation of simple visual assessment (based on leaf/crown transparency [14]) and multispectral indices, eddy covariance techniques, analyses of tree rings, and measurements of foliage/whole-tree photosynthesis rates, as well as sap-flow rates ( [15], and references therein). Those methods may provide detailed information about analyzed dying/healthy trees and, in particular, relevant information on the causal and functional factors driving the mortality of trees in response to intense dry events. However, their application is manpower-, money-, and time-consuming. Moreover, they have clear space-time limitations as they may be applied on a local scale and limited periods.
In a complementary way to field methods, remote sensing techniques may be particularly useful to study climate change-induced phenomena, looking at their effects at a larger spatial and temporal scale. A recent review [16] reports many research works focused on the assessment and monitoring of forest health by using remote sensing techniques from UAV, airborne, and satellite platforms, with passive and active sensor types, multispectral and hyper-spectral technologies, etc. As compared with multispectral information from sensors on manned or unmanned platforms, satellite data allow us to monitor vegetation, computing multispectral vegetation indices which are related to the status of vegetation [17], on a larger spatial and temporal scale. In addition, the use of long time series of satellite data has proved to be suitable to follow space-time dynamics of vegetation and monitor forest decline/mortality in different environments (e.g., birch forest in China [18], oak populations in Iran [19], pine and oak trees in the Mediterranean area [20]).
In the specific context of drought-induced impact on forests, different satellite sensors have been used to monitor and evaluate drought effects on a large scale ( [16,21], and references therein). Among them, instruments aboard the Landsat series are particularly useful for that purpose e.g., [22][23][24] thanks to the right trade-off between temporal revisit time and spatial resolution, as well as their multi-decennial lifetime. On the other hand, the use of suitable satellite data is not enough to extract useful information. Signal natural variability due to phenological phases and weather conditions can mask the effects of climate-change-related stress. Moreover, since site-specific conditions may affect forests on the local scale, analyses should be performed at different spatial levels. Nevertheless, analyses are mainly carried out on a local level or, however, uniquely on the local, regional, or continental scale ( [16], and references therein).
In this work, in order to identify and monitor forests having undergone climate-changerelated stress, a long-time series of Landsat data were exploited to analyze oak stands in Southern Italy both on the local scale (on selected decaying and non-decaying small forests) and on a larger scale (Basilicata region). Two approaches, based on the well-known Normalized Difference Vegetation Index (NDVI, [25]) were used: one based on NDVI values, averaged over the selected forests; another based on the Robust Satellite Techniques (RST, [26,27]). All procedures were implemented on Google Earth Engine (GEE, [28]) since this computing platform has shown to be extraordinarily useful for analyzing long time series of satellite data, particularly in the context of vegetation monitoring (e.g., [18,29,30]).

Material and Methods
In this section, we present data and methodologies we used to evaluate the effects of drought events, occurring around the year 2000, on oak forest areas in Southern Italy. We employed the commonly used Normalized Difference Vegetation Index (NDVI, [25]), based on a long time series of Landsat 5 images (1984Landsat 5 images ( -2011. To this aim, we first analyzed the trends of the NDVI values averaged over selected decaying and similar (e.g., in terms of soil and vegetation characteristics) close non-decaying forests mainly to show what an NDVI-based analysis can highlight. Then, a more refined methodology such as the Robust Satellite Techniques (RST, [26,27]) approach was mainly used to identify in the space-time domain those areas showing anomalous decreases in the NDVI values. Indeed, in different fields of application (e.g., [31][32][33][34][35][36][37]), this robust change detection methodology has proved to be able to identify actual anomalous signals thanks to a preliminary characterization of satellite signals under unperturbed conditions.

Study Area
The Basilicata region (Southern Italy) is characterized by hot dry summers and severe dry winters, where the continental element dominates, as well as humid weather in all seasons, in the west coastal part, due to the Mediterranean Sea influence [38]. Basilicata climate is strongly influenced by orography, with areas characterized by high water deficit up to zones with high precipitation (even higher than 2000 mm/year). On the basis of the Standardized Precipitation Index (SPI, [39]) and the De Martonne aridity index [40] over the period 1993-2004, three areas have been identified [38]: (i) from moderate to extremely humid (Tyrrhenian side), (ii) normal (mainly corresponding to the Basilicata central area), and (iii) dry (Ionian coast and Basilicata Eastern side). The Standardized Precipitation-Evapotranspiration Index (SPEI, [41,42]) over the Basilicata region ( Figure 1) highlights several short dry periods (e.g., in spring 1989, autumn 1990), but it shows a severe and prolonged drought period in the early 2000s. computing platform has shown to be extraordinarily useful for analyzing long time series of satellite data, particularly in the context of vegetation monitoring (e.g., [18,29,30]).

Material and Methods
In this section, we present data and methodologies we used to evaluate the effects of drought events, occurring around the year 2000, on oak forest areas in Southern Italy. We employed the commonly used Normalized Difference Vegetation Index (NDVI, [25]), based on a long time series of Landsat 5 images (1984Landsat 5 images ( -2011. To this aim, we first analyzed the trends of the NDVI values averaged over selected decaying and similar (e.g., in terms of soil and vegetation characteristics) close non-decaying forests mainly to show what an NDVI-based analysis can highlight. Then, a more refined methodology such as the Robust Satellite Techniques (RST, [26,27]) approach was mainly used to identify in the space-time domain those areas showing anomalous decreases in the NDVI values. Indeed, in different fields of application (e.g., [31][32][33][34][35][36][37]), this robust change detection methodology has proved to be able to identify actual anomalous signals thanks to a preliminary characterization of satellite signals under unperturbed conditions.

Study Area
The Basilicata region (Southern Italy) is characterized by hot dry summers and severe dry winters, where the continental element dominates, as well as humid weather in all seasons, in the west coastal part, due to the Mediterranean Sea influence [38]. Basilicata climate is strongly influenced by orography, with areas characterized by high water deficit up to zones with high precipitation (even higher than 2000 mm/year). On the basis of the Standardized Precipitation Index (SPI, [39]) and the De Martonne aridity index [40] over the period 1993-2004, three areas have been identified [38]: (i) from moderate to extremely humid (Tyrrhenian side), (ii) normal (mainly corresponding to the Basilicata central area), and (iii) dry (Ionian coast and Basilicata Eastern side). The Standardized Precipitation-Evapotranspiration Index (SPEI, [41,42]) over the Basilicata region ( Figure  1) highlights several short dry periods (e.g., in spring 1989, autumn 1990), but it shows a severe and prolonged drought period in the early 2000s.  [41]. SPEI categories from [20].
In Basilicata, oak forests represent about 52% of the whole forest area and are present with different structural and composition typologies, from 400-500 to 1200 m a.s.l. [43]. In such a region, field surveys in specific places (e.g., Gorgoglione and San Paolo Albanese, Figures 2-4) have pointed out the decline of oak forests after drought events [9,10,13,20,44,45].       Oak decline in the municipality of Gorgoglione has been evident from visual assessments only since 2003, when about 20% of the forest trees showed dryness of apical branches, bark detachment, and root rot disease [13]. From 2002 to 2004, oak mortality in this area increased from 5% to 10% [10]. Field surveys have been recently carried out in such a mixed forest, which is mainly characterized by Quercus cerris, Quercus pubescens, and other broadleaf species [10]. Field and laboratory analyses of decaying and non-decaying couples of trees [10,13] showed a drought-induced decline since the early 2000s [10].
A similar survey has been carried out in the municipality of San Paolo Albanese in a monospecific forest characterized by Quercus frainetto. Like the Gorgoglione site, the San Paolo Albanese forest has shown shoot dieback, leaf loss and withering, growth decline, and high mortality since the early 2000s [44]. In particular, according to the same authors [44], the onset of dieback was in 2002.
In order to evaluate the stress effects due to the 2000s drought events and the possible climatological forcing, the satellite-based analysis presented in this work was carried out at two different spatial scales: • on a local scale, i.e., analyzing data concerning decaying (D) forests in the municipalities of Gorgoglione and San Paolo Albanese. By comparison, similar (e.g., in terms of soil and vegetation characteristics) not-so-far non-decaying (ND) forests were considered in the municipalities of Accettura-Pietrapertosa and San Paolo Albanese, respectively ( Figures 3 and 4). • on a larger scale, considering oak forests in the Basilicata region.

Satellite Data
To evaluate oak forest changes induced by the series of drought events that occurred around the year 2000, we used data from Thematic Mapper (TM) aboard the Landsat-5 platform. Its spatial resolution (30 m), 16-day revisit period, very long life (active from March 1984 to November 2011), and timespan straggling the period of drought events under study made TM-Landsat 5 the most suitable choice. No data from sensors aboard other Landsat platforms were used in order to avoid the introduction of possible errors due to cross-calibration, which is required by the different spectral response profiles [47]. We considered, as a dataset τ, all available TM-Landsat-5 (Path: 188, Row: 032) images from 1 June to 30 September (i.e., the period including the maximum oak crown extent), from 1984 to 2011. No 2007 cloud-free image is available for the summer season.
In the GEE procedures (described in Section 2.3.2), all available Level 1-Terrain Corrected images, in terms of Top-Of-Atmosphere (TOA) reflectance, belonging to the dataset τ, were considered, for a total amount of 150 images.

Land Cover
The distribution of oak forests in the Basilicata region was obtained from land cover information available on the Italian Institute for Environmental Protection and Research (ISPRA-Istituto Superiore per la Protezione e la Ricerca Ambientale) portal [48]. We downloaded Level 4 of three CORINE Land Cover (CLC, [49]) maps, related to 2000, 2006, and 2012, to take into account all the updated versions within the long-term period investigated. For each CLC map, we selected the informative layers of interest, i.e., 3.1.1.1 (Mediterranean evergreen oak forests) and 3.1.1.2 (Deciduous oak forests). In particular, we pointed at identifying only the areas that maintained the same vegetation cover type in the analyzed period. To this aim, we performed a vector intersection between the three CLC maps to select only the informative layers of interest. In this way, we obtained a final informative layer of (permanent in the considered period) oak forest areas over the Basilicata region ( Figure 5).

Global Forest Change Data
The Hansen [50] Global Forest Change (GFC) dataset was used for comparison with the results achieved by our analysis at the scale of the Basilicata region. It is available in GEE as Hansen Global Forest Change v1.9 product and represents forest change, at 30 m resolution, between 2000 and 2021. The product is based on data from the Enhanced Thematic Mapper Plus (ETM+) sensor aboard Landsat-7 and consists of different bands such as the aggregated forest loss, gain, and year of occurred loss. Information from year of occurred forest loss were used in this work.

NDVI Spatial Average and RST-Based Approach
The analyses carried out in this work are based on the Normalized Difference Vegetation Index (NDVI, [25]). NDVI is a widely used index since it is a powerful indicator of greenness, health, and growth of vegetation and a suitable tool to improve impact assessments of disturbances such as drought [51]. NDVI is computed on the basis of near-infrared (NIR) and visible (VIS) spectral radiances: In this paper, NIR is represented by the TOA reflectance measured in the band B4 (0.76-0.90 μm) of TM-Landsat 5 and VIS is the TOA reflectance measured in the band B3 (0.63-0.69 μm). Three preliminary steps were carried out: (a) the removal of satellite images with georeferencing errors from the dataset τ. The new dataset (τ′) was so made of 148 images; (b) the exclusion of the pixels affected by clouds or cloud shadows from further processing phases; (c) the computation of the NDVI(x,y,t) values, for each pixel centered at (x,y) and not excluded during step b), for each available satellite image acquired at the time t (with t ∈ τ′).
After the three above-mentioned steps, two approaches were used to evaluate the effects of drought events on the oak population.
The first approach we used in this work is based on the computation of the NDVI spatial average for each forest test site. In particular, we calculated mean values of NDVI

Global Forest Change Data
The Hansen [50] Global Forest Change (GFC) dataset was used for comparison with the results achieved by our analysis at the scale of the Basilicata region. It is available in GEE as Hansen Global Forest Change v1.9 product and represents forest change, at 30 m resolution, between 2000 and 2021. The product is based on data from the Enhanced Thematic Mapper Plus (ETM+) sensor aboard Landsat-7 and consists of different bands such as the aggregated forest loss, gain, and year of occurred loss. Information from year of occurred forest loss were used in this work.

NDVI Spatial Average and RST-Based Approach
The analyses carried out in this work are based on the Normalized Difference Vegetation Index (NDVI, [25]). NDVI is a widely used index since it is a powerful indicator of greenness, health, and growth of vegetation and a suitable tool to improve impact assessments of disturbances such as drought [51]. NDVI is computed on the basis of near-infrared (NIR) and visible (VIS) spectral radiances: In this paper, NIR is represented by the TOA reflectance measured in the band B4 (0.76-0.90 µm) of TM-Landsat 5 and VIS is the TOA reflectance measured in the band B3 (0.63-0.69 µm). Three preliminary steps were carried out: (a) the removal of satellite images with georeferencing errors from the dataset τ. The new dataset (τ ) was so made of 148 images; (b) the exclusion of the pixels affected by clouds or cloud shadows from further processing phases; (c) the computation of the NDVI(x,y,t) values, for each pixel centered at (x,y) and not excluded during step b), for each available satellite image acquired at the time t (with t ∈ τ ).
After the three above-mentioned steps, two approaches were used to evaluate the effects of drought events on the oak population.
The first approach we used in this work is based on the computation of the NDVI spatial average for each forest test site. In particular, we calculated mean values of NDVI (x,y,t) within the D areas of Gorgoglione ("GORG") and San Paolo Albanese ("SPA") as well as the ND forests of Accettura-Pietrapertosa ("AP") and San Paolo Albanese (Figures 3 and 4), in order to generate NDVI GORG_D (t), NDVI SPA_D (t), NDVI AP_ND (t), and NDVI SPA_ND (t), respectively, for each available satellite image acquired at the time t (t ∈ τ ).
The second approach we used is based on the Robust Satellite Techniques (RST, [26,27]). Such a methodology has been successfully applied both to polar and geostationary satellite sensor data in different fields of application such as volcanic (e.g., [52]), seismic (e.g., [53]), hydrological [54,55], and fire [37] risks, as well as in the context of cloudy-radiance detection [56], Saharan dust (e.g., [57]), gas flaring (e.g., [33]), and sea quality (e.g., [34,58]) monitoring. RST is a change detection methodology that identifies as "anomalous" the signal that deviates significantly from its behavior in the normal conditions of the specific place and time of observation. It relies on the Absolutely Local Index of Change of the Environment (ALICE), so defined: where V(x,y,t) is the value of the V signal, measured at the place with coordinates (x,y) and at the time t. µ V (x,y) and σ V (x,y) are the temporal mean and standard deviation, respectively, and represent the expected value and the normal variability of the signal V(x,y,t) itself. We will refer to them as "reference fields" henceforth. They are computed on the basis of a multi-year dataset of cloud-free satellite records, which are collected under the same observational conditions (e.g., same platform, time of day, month/season of the year) as the image to be processed. By construction, the index in Equation (2) • reduces the effects of known sources of natural (e.g., topography, land cover) and observational (e.g., solar illumination, satellite view angle) noise; • takes account of the effects of the residual noise due to the non-predictable signal variability in the space-time domain (e.g., changes in soil and atmospheric moisture); • is intrinsically robust against false alarm proliferation.
Moreover, it should be highlighted that the index (2) is a standardized variable, for its inherent design, that should show a Gaussian behavior (i.e., mean~0 and standard deviation~1). The higher the absolute value measured, the lower the probability of occurrence. In particular, the probability of occurrence of an anomalous signal decreases from about 4.5% for |⊗ V | > 2 (level of statistical significance) to 0.3% for |⊗ V | > 3 and even lower values for increased absolute values of ⊗ V [59].
In the RST-based analysis we carried out, V(x,y,t) is not represented by the simple NDVI(x,y,t) value of each pixel centered at (x,y) on the TM-Landsat-5 image acquired at the time t (with t ∈ τ ). In order to better characterize the signal corresponding to the maximum crown extent as well as reduce cloud contamination (possibly missed in the cloud detection phase) and other effects due to different observational conditions (e.g., atmospheric transmittance, satellite view angle), the Maximum Value Composite (MVC, [60]) of NDVI was indeed considered as V(x,y,t). Therefore, this second approach may be summarized in the following phases: • the computation of NDVI MVC (x,y,t ) for each pixel (x,y) and each year t (t ∈ τ ), retaining the highest NDVI(x,y,t) value on a 4-month basis (June, July, August, and September); • the calculation of reference fields, µ NDVI_MVC (x,y) and σ NDVI_MVC (x,y), considering all NDVI MVC (x,y,t ) with t ∈ [1984][1985][1986][1987][1988][1989][1990][1991][1992][1993][1994][1995][1996][1997] (i.e., considered as an unperturbed or less perturbed period); • the computation of the ALICE index (⊗ NDVI_MVC (x,y,t )) for each pixel (x,y) and t' ∈ , following Equation (3): ⊗ NDVI_MVC (x,y,t ) gives the local excess of the NDVI MVC (x,y,t ) signal compared with its historical mean value, µ NDVI_MVC (x,y), and weighted by its historical variability σ NDVI_MVC (x,y), at the considered location (x,y) and time of observation t . In the case of decaying forests, anomalous decreasing values of NDVI MVC (x,y,t ) and, consequently, negative ⊗ NDVI_MVC (x,y,t ) are expected. In particular, the more anomalous NDVI MVC (x,y,t ), the more negative ⊗ NDVI_MVC (x,y,t ).
In our analysis, we wanted to highlight phenomena with a low probability of occurrence and, for this reason, taking into consideration what is said before, we considered as "strongly anomalous" those pixels with ⊗ NDVI_MVC (x,y,t ) ≤ −3.

Procedure Implementation on Google Earth Engine
Landsat-TM Data Pre-Processing, NDVI Trend Extraction, and RST-Based Processing The two approaches before described were implemented on the GEE platform [28]. TOA reflectances of Landsat 5 TM Collection 1 Tier 2 were used in our investigation. As already mentioned, we analyzed the scenes collected over the Basilicata region (Path: 188, Row: 032), in June, July, August, and September, from 1984 to 2011. The initial 150 available scenes were reduced to 148 after applying the filter requiring a high geometric accuracy (i.e., "GEOMETRIC_RMSE_VERIFY" < 1).
To filter out the pixels contaminated by clouds, the TM-Landsat quality band was used to select only clear pixels (i.e., "BQA" = 672). Moreover, aimed at removing residual undetected reflecting clouds, the RST-based OCA approach (One-channel Cloud-detection Approach, [56]) was used by exploiting channel B1 (0.45-0.52 µm). Reference fields, µ B1 (x,y) and σ B1 (x,y), were computed on a monthly basis and a multi-year dataset , following RST rules. Then, we applied Equation (2) with V(x,y,t) ≡ B1(x,y,t): Pixels with ⊗ B1 (x,y,t) ≥ 1 were considered contaminated by clouds. Finally, the locations affected by cloud shadows, which were undetected by the CFMask (C Function of Mask, [61]) algorithm, were identified by taking into account the information on Sun azimuth angle stored in the metadata (i.e., "SUN_AZIMUTH") file for each TM image.
For each pixel that was not affected by clouds/shadows, the NDVI index (Equation (1)) was computed by using band B4 (0.76-0.90 µm) and band B3 (0.63-0.69 µm). After the ingestion of the vector layers of D and ND forests, we computed the NDVI spatial average for all test sites and then plotted them as a time series (see next section, Figure 6). Finally, as above-mentioned, we implemented the RST approach on the MVC of NDVI, rather than the simple NDVI. So, for each year t (t ∈ ) and clear pixel, the highest value of the NDVI recorded from June to September was used to generate a synthetic image that was representative of the annual MVC, i.e., NDVI MVC (x,y,t ). A total of 27 NDVI MVC images were generated, from 1984 to 2011, being 2007 not available. Using NDVI MVC images from 1984 to 1997, the RST reference fields, µ NDVI_MVC (x,y) and σ NDVI_MVC (x,y), were computed. In the end, for each pixel of the annual synthetic NDVI MVC images, the ALICE index ⊗ NDVI_MVC (x,y,t ) was computed as described in Equation (3). After the ingestion of the vector layer of Basilicata oak forests as well, the following RST-based products were generated: Different drought-tolerance of species should be considered among the possible reasons of the trend differences observed between the analyzed oak forests. In fact, it may be supposed a different drought response of a mixed stand, such as the Accettura-Pietrapertosa/Gorgoglione forests, respect to a monospecific stand such as in San Paolo Albanese.
It should be highlighted that the above-described trends of the NDVI spatial average over D/ND forests in Gorgoglione/Accettura-Pietrapertosa and San Paolo Albanese are in agreement with findings obtained in a recent study [20] over the same areas using the maximum annual NDVI value of MODIS aggregated products, over the period 2001-2019.
This first analysis allowed us to highlight the raising trends of the forests on the whole as well as the temporary decrease around the year 2000.
Concerning the RST-based approach, to highlight vegetation areas with a worsening of health over the analyzed period, strongly anomalous NDVIMVC(x,y,t′) values (⨂NDVI_MVC (x,y,t′) ≤ −3) were considered. As described in the implementation section, such anomalous pixels were counted per year and over each investigated area (D and ND forests in Gorgoglione/Accettura-Pietrapertosa and San Paolo Albanese as well as oak forests in the Basilicata region). The results are shown in Figure 7a (Gorgoglione/Accettura-Pietrapertosa), Figure 7b (San Paolo Albanese), and Figure 9 (Basilicata). They are presented as percentages since areas have a different extension and, in order to make even small variations clear, the scale of each graphic is adapted to the data of each test area.
The RST-based analysis on the local scale ( Figure 7) shows that: as above-mentioned ( Figure 6), a different initial health status of forests seems to exist even before 2000s. This condition, together with local factors, could be responsible For each year and within the Basilicata oak forest informative layer, we computed the annual forest trend, by counting the pixels with value 1 in the "lossyear" band.

Results and Discussion
A comparative analysis between D and ND forests was made considering the results of the first approach, based on the NDVI spatial average, in order to highlight possible different trends, particularly in correspondence with drought events in the 2000s. Figure 6a reports the values of NDVI GORG_D and NDVI AP_ND for each available TM-Landsat-5 summer image, from 1984 to 2011. It is clear that since the 1980s, the ND forest has had greater values than the D one. This indicates that the ND forest has always been more vigorous and/or with dense canopy.
Their variations seem to be synchronous, indicating a common forcing factor, and, in fact, trend lines run almost in parallel. The NDVI values of the two forests do not show any clear evidence of a convergent or divergent trend. Both of them record a decrease in the early 2000s, but the time series shows minimum values in other periods too (NDVI GORG_D is 0.466 on 9 September 1990; NDVI GORG_ND is 0.539 on 20 August 2006). On the whole, the NDVI values show a rising trend and, after their abrupt decrease around 2000, they try to reach again, in the following years, the values of the last 1990s. Although, as said before, vegetation response to dryness periods may be different in time and space, it is worth noting the likeness between such time series and SPEI-12 ( Figure 1). In particular, some plot elements (e.g., the prolonged signal depression around the year 2000 and the subsequent increase, followed by a new signal decrease) are similar. Figure 6b reports the trends for NDVI SPA_D and NDVI SPA_ND . Likewise, Accettura-Pietrapertosa/Gorgoglione test sites, the two time-series of NDVI related to San Paolo Albanese forests seem to be synchronous, but on the whole, the trend lines show a quite Different drought-tolerance of species should be considered among the possible reasons of the trend differences observed between the analyzed oak forests. In fact, it may be supposed a different drought response of a mixed stand, such as the Accettura-Pietrapertosa/Gorgoglione forests, respect to a monospecific stand such as in San Paolo Albanese.
It should be highlighted that the above-described trends of the NDVI spatial average over D/ND forests in Gorgoglione/Accettura-Pietrapertosa and San Paolo Albanese are in agreement with findings obtained in a recent study [20] over the same areas using the maximum annual NDVI value of MODIS aggregated products, over the period 2001-2019.
This first analysis allowed us to highlight the raising trends of the forests on the whole as well as the temporary decrease around the year 2000.
Concerning the RST-based approach, to highlight vegetation areas with a worsening of health over the analyzed period, strongly anomalous NDVI MVC (x,y,t ) values (⊗ NDVI_MVC (x,y,t ) ≤ −3) were considered. As described in the implementation section, such anomalous pixels were counted per year and over each investigated area (D and ND forests in Gorgoglione/Accettura-Pietrapertosa and San Paolo Albanese as well as oak forests in the Basilicata region). The results are shown in Figure 7a   On the whole, we may suppose that the observed peaks in the early 2000s could be a first direct consequence of loss in tree vigor and mortality due to a compromised hydraulic functionality. Some authors [13] report lumen vessels filled with tyloses in 1999, 2000, and 2001 in some analyzed trees of the Gorgoglione forest, indicating a decreased xylem conductivity and an increased hydraulic vulnerability, i.e., a condition contributing to tree On the whole, we may suppose that the observed peaks in the early 2000s could be a first direct consequence of loss in tree vigor and mortality due to a compromised hydraulic functionality. Some authors [13] report lumen vessels filled with tyloses in 1999, 2000, and 2001 in some analyzed trees of the Gorgoglione forest, indicating a decreased xylem conductivity and an increased hydraulic vulnerability, i.e., a condition contributing to tree vigor reduction and death [62]. The same forest was specifically investigated by the same authors [13] in 2013 also from a phytopathological point of view. No biotic factors were identified as the triggering cause of the forest dieback. In fact, the parasitic attacks observed are mostly attributable to pathogens of weakness, which settle in plants that are generally already suffering [13]. This is also in agreement with more general observations, which usually recognize the pathogens associated with the oak decline in the Mediterranean area as secondary biotic factors [9].
Field sampling and laboratory analyses in San Paolo Albanese D forest highlighted growth decline of the stand in response to warm and dry spring conditions with the onset of dieback in 2002 [44]. Declining trees displayed a sharp growth drop after this onset with reductions of 49% [44].
The analysis of the ALICE ⊗ NDVI_MVC (x,y,t ) maps, for each year t ∈ [1984-2011], allowed us to perform a spatial investigation. In particular, the worsening of oak forest health status is evident around the year 2000 on a local scale comparing 1999, 2000, and 2001 maps (Figure 8). The differences are in terms of both extension and anomaly intensity (i.e., lower ALICE values). The same analysis allowed us to highlight also that the Accettura-Pietrapertosa forest is "non-decaying" looking at it on the whole, but it has small decaying parts inside (Figure 8b).
The analysis of the ALICE ⨂NDVI_MVC (x,y,t′) maps, for each year t′ ∈ [1984-2011], allowed us to perform a spatial investigation. In particular, the worsening of oak forest health status is evident around the year 2000 on a local scale comparing 1999, 2000, and 2001 maps (Figure 8). The differences are in terms of both extension and anomaly intensity (i.e., lower ALICE values). The same analysis allowed us to highlight also that the Accettura-Pietrapertosa forest is "non-decaying" looking at it on the whole, but it has small decaying parts inside (Figure 8b).  For comparison, Figure 9 also shows the annual trend of the GFC loss product [50] from 2001 to 2011. It is worth noting that the GFC product and RST results partially differentiate since the Hansen product reports vegetation loss, while the RST-based analysis highlights statistically significant ("anomalous") variations of vegetation. This dissimilarity is reflected in the different percentages reported by RST and GFC in the plot of Figure  9. However, both of them measure parameters affected by the same climatological forcing, which is responsible for drought and forest decaying, i.e., favorable factors for forest fire For comparison, Figure 9 also shows the annual trend of the GFC loss product [50] from 2001 to 2011. It is worth noting that the GFC product and RST results partially differentiate since the Hansen product reports vegetation loss, while the RST-based analysis highlights statistically significant ("anomalous") variations of vegetation. This dissimilarity is reflected in the different percentages reported by RST and GFC in the plot of Figure 9. However, both of them measure parameters affected by the same climatological forcing, which is responsible for drought and forest decaying, i.e., favorable factors for forest fire increase as well as selection of tree cutting. Such a common element is confirmed by the good correlation between the two trends (0.67 Pearson correlation coefficient) and by the two important peaks (recorded in 2001 and 2004) in both trends.
The worsening of oak forest health status is evident around the year 2000 also at the Basilicata level, through the analysis of the RST maps ( Figure 10). For comparison, Figure 9 also shows the annual trend of the GFC loss product [50] from 2001 to 2011. It is worth noting that the GFC product and RST results partially differentiate since the Hansen product reports vegetation loss, while the RST-based analysis highlights statistically significant ("anomalous") variations of vegetation. This dissimilarity is reflected in the different percentages reported by RST and GFC in the plot of Figure  9. However, both of them measure parameters affected by the same climatological forcing, which is responsible for drought and forest decaying, i.e., favorable factors for forest fire increase as well as selection of tree cutting. Such a common element is confirmed by the good correlation between the two trends (0.67 Pearson correlation coefficient) and by the two important peaks (recorded in 2001 and 2004) in both trends.
The worsening of oak forest health status is evident around the year 2000 also at the Basilicata level, through the analysis of the RST maps ( Figure 10).  On the other hand, the map reporting pixel-by-pixel the count of anomalous NDVI MVC (x,y,t ) values (⊗ NDVI_MVC (x,y,t ) ≤ −3), from 1984 to 2011, permitted us to highlight other decaying oak forests in Basilicata. Figure 11 shows an example of decaying areas in the municipality of Garaguso, near the Salandrella torrent and internal paths, in agreement with the forest report [63]. It is worth underlining that the patchy aspect of such maps reflects the possible interference of local biotic (e.g., forest structure and composition) and abiotic (e.g., topography and soil) factors in the drought effects on forest changes [22].
NDVIMVC(x,y,t′) values (⨂NDVI_MVC(x,y,t′) ≤ −3), from 1984 to 2011, permitted us to highlight other decaying oak forests in Basilicata. Figure 11 shows an example of decaying areas in the municipality of Garaguso, near the Salandrella torrent and internal paths, in agreement with the forest report [63]. It is worth underlining that the patchy aspect of such maps reflects the possible interference of local biotic (e.g., forest structure and composition) and abiotic (e.g., topography and soil) factors in the drought effects on forest changes [22].

Conclusions
Two approaches were used to analyze possible climate change-related stress over oak forests in the Basilicata region. The analysis based on the NDVI spatial average on D and similar close ND forests gave us information about the initial different health status of the test areas as well as the common rising NDVI trend from 1994 to 2011 despite a clear multi-year decrease in values around the year 2000.
Differently from the first approach, the RST-based analysis takes account of signal natural variability by its inherent design and was able to better highlight differences of the forests both in space and time domains. The analysis performed over time at different spatial levels allowed us to observe stress effects on a local and regional scale. Site-specific peaks of the number of strongly anomalous (⨂NDVI_MVC(x,y,t′) ≤ −3) pixels could be the result of adverse local effects, while common peaks (e.g., in 2001) are expected in the case of a climatological forcing. In particular, the RST-based analysis permitted us to emphasize different stress impacts in the test case areas. With similar trees and under the same stressed conditions, the initial worse health status of the vegetation in the D forest probably caused a more marked effect that was pointed out by a greater percentage of strongly anomalous pixels. At the regional level, different oak species and local conditions may have contributed to a different response to drought. A quite important result is the evidence of strong anomalous values of NDVIMVC already in 2001 in the Gorgoglione forest. In fact, the combination of suitable satellite data and a sensitive methodology allowed us to clearly highlight a decaying status of vegetation two years before the first visual field evidence [13].
On the other hand, the spatial analysis was able to highlight variations in extension and anomaly intensity over time. Moreover, it identified decaying areas (even within an overall non-decaying forest) to be focused on. Such affected areas could be studied, for example, with field-based techniques. In fact, a multi-level approach, from lab/field experiments to long time series of satellite data analyses, should be preferred in order to

Conclusions
Two approaches were used to analyze possible climate change-related stress over oak forests in the Basilicata region. The analysis based on the NDVI spatial average on D and similar close ND forests gave us information about the initial different health status of the test areas as well as the common rising NDVI trend from 1994 to 2011 despite a clear multi-year decrease in values around the year 2000.
Differently from the first approach, the RST-based analysis takes account of signal natural variability by its inherent design and was able to better highlight differences of the forests both in space and time domains. The analysis performed over time at different spatial levels allowed us to observe stress effects on a local and regional scale. Site-specific peaks of the number of strongly anomalous (⊗ NDVI_MVC (x,y,t ) ≤ −3) pixels could be the result of adverse local effects, while common peaks (e.g., in 2001) are expected in the case of a climatological forcing. In particular, the RST-based analysis permitted us to emphasize different stress impacts in the test case areas. With similar trees and under the same stressed conditions, the initial worse health status of the vegetation in the D forest probably caused a more marked effect that was pointed out by a greater percentage of strongly anomalous pixels. At the regional level, different oak species and local conditions may have contributed to a different response to drought. A quite important result is the evidence of strong anomalous values of NDVI MVC already in 2001 in the Gorgoglione forest. In fact, the combination of suitable satellite data and a sensitive methodology allowed us to clearly highlight a decaying status of vegetation two years before the first visual field evidence [13].
On the other hand, the spatial analysis was able to highlight variations in extension and anomaly intensity over time. Moreover, it identified decaying areas (even within an overall non-decaying forest) to be focused on. Such affected areas could be studied, for example, with field-based techniques. In fact, a multi-level approach, from lab/field experiments to long time series of satellite data analyses, should be preferred in order to highlight possible trends as well as the actual space-time scale of events, as suggested by [64,65]. In particular, some authors [20] suggest that the combination of analyses based on remote sensing and tree-ring data may be a useful tool to better understand droughtinduced phenomena in forests. Indeed, this combined approach could be really useful in the case of time lag between the remote sensing signal and the climatic stresses due to, for example, the likely peculiar site-specific conditions such as low tree density or non-uniform canopy coverage.
Finally, it should be underlined the peculiarity of this work, which is based on the combination of three elements: • suitable satellite data, in terms of spatial resolution, temporal repetition, and lifetime span; • the RST-based analysis, which is fully exportable to different sensors (as already proved in different contexts) and independent on the specific forest type thanks to its pixel-based signal characterization (i.e., reference fields); • the potential of GEE, which allows us to analyze long time series of satellite data without computational efforts, from a local to a larger scale.
The combination of these three elements, here used in the case of the Basilicata region, may be exploited, in general, to monitor long-term forest health, from the local to global scale, as well as to timely detect climate change forcing. In particular, the identification of more drought-vulnerable areas over large areas could be useful to forest managers and policymakers [66], while monitoring of annual forest decline could help ecologists to evaluate climate change impact on vegetation [22]. Finally, it is worth noting that the RST exportability allows us to employ other satellite sensors. Although Landsat satellite is still the most employed mission in this kind of study [16], RST could be easily applied to Sentinel-2 imagery, which has greater repetition time (5 days) and higher spatial resolution (up to 10 m). This improvement could help, for example, in better integrating field and remote sensing methodologies since estimates of forest variables depend on how field data can be linked to satellite images [67].