Recurrence Analysis of Vegetation Indices for Highlighting the Ecosystem Response to Drought Events: An Application to the Amazon Forest

Forests are important in sequestering CO2 and therefore play a significant role in climate change. However, the CO2 cycle is conditioned by drought events that alter the rate of photosynthesis, which is the principal physiological action of plants in transforming CO2 into biological energy. This study applied recurrence quantification analysis (RQA) to describe the evolution of photosynthesis-related indices to highlight disturbance alterations produced by the Atlantic Multidecadal Oscillation (AMO, years 2005 and 2010) and the El Niño-Southern Oscillation (ENSO, year 2015) in the Amazon forest. The analysis was carried out using Moderate Resolution Imaging Spectroradiometer (MODIS) images to build time series of the enhanced vegetation index (EVI), the normalized difference water index (NDWI), and the land surface temperature (LST) covering the period 2001–2018. The results did not show significant variations produced by AMO throughout the study area, while a disruption due to the global warming phase linked to the extreme ENSO event occurred, and the forest was able to recover. In addition, spatial differences in the response of the forest to the ENSO event were found. These findings show that the application of RQA to the time series of vegetation indices supports the evaluation of the forest ecosystem response to disruptive events. This approach provides information on the capacity of the forest to recover after a disruptive event and, therefore is useful to estimate the resilience of this particular ecosystem.


Introduction
Terrestrial ecosystems such as forests absorb approximately 60 Gt of carbon annually [1]. Simultaneously, autotrophic and heterotrophic organisms release approximately the same amount of carbon back into the atmosphere, thus closing the terrestrial carbon cycle. Small alterations in the terrestrial carbon budget are likely to have a significant impact on atmospheric CO 2 concentrations [2,3]. There has been substantial warming in the northern high latitudes, which has been associated with changes in the inter-annual variability of the global carbon cycle [4,5]. Documenting and interpreting changes in the ecosystem functions that influence the carbon cycle globally is critical for understanding the associations and feedback between temperature, land cover, and atmospheric CO 2 [6,7]. and the land surface temperature (LST). These are considered as a proxy of the major vegetation variables such as chlorophyll, leaf water content, and temperature, which influence photosynthesis. The analysis of the synchronization of the vegetation indices is expected to effectively identify the types of response of the forest ecosystem to disturbances produced by extreme climate events, highlighting changes in the recurrent pattern of periodic and irregular cyclicities that characterize the ecological function.

The Study Area
The study area is a large area of the Amazon forest. The Amazon forest is a tropical ecosystem characterized by wet conditions; however, it has experienced recurrent and large-scale droughts over the last few years [44,45].  [44,46,47].

Data Collection and Vegetation Indices
With the aim of analyzing vegetation indices in the study area, data were collected from Moderate Resolution Imaging Spectroradiometer (MODIS) images provided by the National Aeronautics and Space Administration (NASA) with a high daily frequency. In particular, the MODIS13A2 product was used to build EVI and NDWI time series profiles. These are composite images, so for each pixel, a value was selected from all the acquisitions within the 16-day composite period. The algorithm for this product chooses the best available pixel value from all the acquisitions of the 16-day period. The criteria used are: (1) low clouds, (2) low view angle, and (3) the highest normalized difference vegetation index (NDVI)/EVI value [48] For both indices, the bandwidth of the near infrared (NIR) was used to amplify the information contained in the RED band (which is sensitive to the content of chloroplasts), and in the band of short wave infrared (SWIR) (which is sensitive to the water content). The BLUE band is used for atmospheric correction in EVI, which is one important feature for studies in the Amazon basin where seasonal burning of pasture and forest takes place throughout the dry season, either for agricultural purposes (land clearing) or during natural fire events. EVI has been used for the study of temperate forests and is much less sensitive to aerosols (from biomass burning) and saturation [49][50][51][52]. In particular, EVI is defined as the normalized ratio between the RED, NIR, and BLUE bands as follows: This index is dimensionless and can assume values between −1 and +1, but for vegetation, generally values between 0 (no vegetation) and 1 occur, depending on the vegetation cover. A low value represents areas with low or no vegetation, implying that the vegetation is in a stressful situation; on the other hand, high or very high values reflect a situation of strong photosynthetic activity and the high presence of biomass. EVI addresses the photosynthetically active chloroplasts, but it is not strongly related to the water content. However, this index alone has some limitations when employed for photosynthesis analysis because it is strongly related to chlorophyll content, and to a lower extent, to water content in the vegetation [53][54][55][56][57]. The water is an important variable affecting photosynthesis because it influences the opening of the leaves stomata that regulate gaseous exchanges. In fact, the water content is a crucial factor for the regulation of both the rate of photosynthesis and the speed of production and the development of new leaves [52]. Therefore, when using vegetation indices like a proxy to study photosynthesis function on vegetation, it should be considered that [10,58]: • each plant species has its own relationship of chlorophyll content and water content; • a decrease in chlorophyll content does not imply a decrease in water content; and • a decrease in volumetric water content does not imply a decrease in chlorophyll content.
To overcome this limitation, the analysis used the NDWI, which is a remote sensing-based index that is sensitive to the change in the water content of leaves [10,54]. NDWI provides information on the impact that stress can have on water content at the leaf level [59,60]. By employing the NDWI, it is possible to integrate the presence of green photosynthetically active biomass to the presence of water, and check that plants did not get hydric stress. It is defined as follows: It ranges between −1 to 1: the higher the water content, the higher the value of NDWI, while NDWI decreases in periods of water stress. Combined with EVI, the NDWI is used to evaluate the water content (and thus the stress) in plants.
The combined recurrent pattern of the EVI and NDWI time series constitutes a good tool for the evaluation of vegetation cover degradation in relation to external pressure factors and the photosynthesis of a certain area. In order to analyze the climate effect, LST profiles using the MOD11A2 product have also been built. Mainly, LST is treated as a direct effect of high temperature and drought events and as an indirect cause of plant stress. This was necessary because we did not have spatially-representative meteorological data. The MOD11A2 Version 6 product provides an average 8-day per-pixel land surface temperature and emissivity (LST&E) with a 1 km spatial resolution. Each pixel value in the MOD11A2 is a simple average of all the corresponding MOD11A1 LST pixels collected within that 8-day period [48]. LST time series were resampled from an eight days frequency to 16 days frequency by averaging two MOD11A2 time windows corresponding to the 16-day time window of MOD13A2. Since MODIS images can produce data of different quality, data with low quality (i.e., pixel not produced due to cloud effects and pixel not produced primarily due to reasons other than cloud) were searched and corrected for missing and erroneous data by reviewing the quality assurance flags that were provided together with the data. Such data were here replaced using interpolation of the first good quality temporally adjacent data at the previous and following time steps.
Three-time profiles covering 16 days of average values of EVI, NDWI, and LST were thus constructed from January 2001 to December 2018, both for undisturbed areas. The areas analyzed were selected by digitalizing homogeneous forested areas by means of ortophoto interpretation excluding areas affected by human activities. This operation was carried out using QGIS software, which allowed the visualization of the Google Earth and Bing Maps orthophoto. Next, we verified that the area selected did not change in the years 2005, 2010, and 2015 using historical images from Google Earth. This step was necessary to exclude the effect of human activities and other pressures on vegetation that may affect the ecological function. The time series were extracted using the Application for Extracting and Exploring Analysis Ready Samples (AppEEARS) provided by the U.S. Geological Survey. Furthermore, meteorological data available from "Wunderground" were analyzed for the period January 2001-December 2018. Specifically, the temperature and the humidity were collected using Santarem Station (IMOJUD2), which is located in Brazil [61]. Time series of temperature and humidity were built using the same frequency of the EVI and NDWI time series (16 days) to allow for a comparison between vegetation indices and meteorological data.

Recurrence Analysis
Recurrence of vegetation states over a season is obvious and their investigation using recurrence analysis seems promising. A recurrence means that the vegetation system recurs to a previous state (up to a small variation). For example, the seasonality effect in the vegetation time series causes periodic recurrences. A change of the seasonality can be analyzed with a recurrence plot analysis that uncovers typical or atypical recurrence patterns [62][63][64]. A recurrence plot (RP) is an advanced technique of nonlinear data analysis. It is a visual representation of all pairs of time points (i, j) at which the dynamical system recurs to a previous state x i : Here, we study the relationship between the EVI, NDWI, and LST time series by extensions of the recurrence plot approach. The cross recurrence plot (CRP) checks the co-occurrence of similar states of the EVI, NDWI and LST, whereas the joint recurrence plot (JRP) analyzes the simultaneous occurrence of recurrences in the time series [33,65]. The two methods compare different aspects of the considered time series.
Mathematically, a CRP between two variables x i and y j is calculated by: and a JRP by the element-wise multiplication of the individual recurrence plots of each system x i and y i : with the states x i ,y i ∈ R representing the time series of two indices coupled as the EVI-NDWI, EVI-LST, NDWI-LST; i,j = 1 . . . N; N is the number of considered states x i ,y i ; is a threshold distance (separate thresholds for signals x and y in the JRP, chosen in a way that ensures the same temporal interval in the evaluation of the each RP); ||· ||s a norm; and θ(·) is the Heaviside function [66].
In the CRP analysis, the comparison is made between similar systems. It checks for the occurrences of similar states (values) in the analyzed signals. The main measure of similarity is the line of synchronization (LOS), which is the analogue to the main diagonal in the RP (the line of identity). If both signals are equal, the LOS is a straight line and is the main diagonal in the plot. Slight deviations in one of the systems cause a distortion of the LOS (e.g., interruptions or bendings, Figure 3c). For signals that are uncorrelated, the LOS completely vanishes. The LOS can be used to identify complete synchronization or time distortions (e.g., two periodic signals of different frequencies would have a LOS with a slope that corresponds to the ratio of the frequencies) [67]. In contrast, the JRP always has a main diagonal, because it is based on the individual RPs of each signal. Here, the comparison was performed by calculating the loss of recurrence points in the JRP relative to the individual RPs. If the considered signals share the identical recurrence patterns, the JRP is the same as the individual RPs (i.e., all recurrence points in the original RPs are preserved in the JRP). The recurrence patterns can even be the same when the signals differ at a first glance (e.g., in the first two moments). This allows for checking for generalized synchronization. In Figure 1 shows an illustration of the cross and joint recurrence plots.
In order to check for complete and generalized synchronization between the EVI, NDWI, and LST time series, we calculated the CRP and JRP between them and quantified them using the following measures. As the thresholds ε x and ε y were chosen to preserve the same time of the RPs for signals x and y, the recurrence rate (RR) for both signals were identical.
For the identification and quantification of complete synchronization, we used the F-recurrence measure given by 1 N i C i,i , which quantifies the density of points along the main diagonal.
The F-recurrence parameter is 1 for complete synchronization and tends to 0 for completely desynchronized systems. This CRP and JRP analyses were combined with a sliding window approach in order to investigate the temporal change of the synchronization. We applied overlapping, sliding temporal windows with a one-year length to the original time series to calculate the F-recurrence measure before the El Niño event (from 2012 to 2014) and after (from 2016 to 2018); and before the two AMO events (for the first event from 2002 to 2004, for the second event from 2007 to 2009) and after (from 2006 to 2008 and from 2011 to 2013, respectively). For the recurrence threshold, we used an empirically defined threshold that ensured a fixed RR = 0.2. The selected value does not qualitatively affect the recurrence quantification analysis (RQA) measures [61,62,67]. The CRP and JRP analyses were carried out in MATLAB using the CRP Toolbox available at http://tocsy.pik-potsdam.de/CRPtoolbox/.
Next, we employed the Fisher test to evaluate the Gaussian distribution of the F-recurrence parameter values as well as the Student's T-test to assess if the means of the F-recurrence parameter distributions before and after the ENSO and AMO events were the same. identity). If both signals are equal, the LOS is a straight line and is the main diagonal in the plot. Slight 210 deviations in one of the systems cause a distortion of the LOS (e.g., interruptions or bendings, Figure   211 3c). For signals that are uncorrelated, the LOS completely vanishes. The LOS can be used to identify 212 complete synchronization or time distortions (e.g., two periodic signals of different frequencies 213 would have a LOS with a slope that corresponds to the ratio of the frequencies) [67]. In contrast, the 214 JRP always has a main diagonal, because it is based on the individual RPs of each signal. Here, the 215 comparison was performed by calculating the loss of recurrence points in the JRP relative to the 216 individual RPs. If the considered signals share the identical recurrence patterns, the JRP is the same 217 as the individual RPs (i.e., all recurrence points in the original RPs are preserved in the JRP). The 218 recurrence patterns can even be the same when the signals differ at a first glance (e.g., in the first two  Figure 2 shows the areas of the Amazon forest selected for the present analysis. These areas showed similar spatial patterns with no changes in structure and composition during the time of the analysis (i.e., from 2001 to 2018). We can therefore rule out that changes in these areas were caused by human activities or fires. The main variation in the vegetation indices can thus be linked to the change in the ecological function of the vegetation forest due to the drought effect.              The profiles give an indication of the temporal variation of forest functions linked with photosynthesis in the analyzed area of the Amazon forest. The figures show that the EVI, NDWI, and LST are characterized by both periodic (represented by the wave recurrent pattern) and non-periodic (represented by peaks characterizing each wave) components.

Time Series
Time series of maximum temperature and humidity were also constructed from 2001 to 2018 with 16 day frequencies (Figure 8) in order to analyze the link between the ecological function of the forest and climate conditions. These profiles help to explore the specific meteorological data measured at the station of the investigated area and provide a way of assessing the meteorological variation in the area. As found for the indices, these time series were characterized by both periodic and non-periodic components.

Synchronization and Evolution Analysis
The individual RPs of the EVI and NDWI showed periodic activities representing the seasonal cycle of the forest ecosystem. The JRP between the EVI and NDWI remained as the periodic recurrence pattern over the years, with a regular and periodic temporal profile without particular variations or disruptions. This confirms the simultaneous recurrent pattern of the vegetation dynamics and photosynthetic activity within the undisturbed area ( Figure 9).   In contrast, the RP of the LST revealed a disruptive event in 2015. Subsequently, the JRP showed periodic recurrences between the EVI and LST over the years, but with disruptions in the years 2015-2016 ( Figure 10). The same situation was present in the JRP between the NDWI and LST. This means that some disturbing event such as an extreme climate event was causing temporal sudden changes in the system that characterizes the LST. In the last period (from 2016 to 2018), the JRP looked very similar to the individual recurrence plots of the first period (from 2001 to 2014) (i.e., the recurrence structure of both signals was the same). Therefore, the system therefore recovered its status afterward, showing a good resilience as the JRP plot maintained a similar structure after the disruption.  These disruptions can be linked to factors that act on short time scales [65]. By analyzing the recurrent pattern of the time series of T max and H avg obtained from the meteorological station, the same disruption could be observed in the same period ( Figure 10). Although these meteorological data are not representative of the entire study area, they may provide greater certainty that the anomaly recorded in the LST for 2015 and 2016 was not stochastic or related to issues in the construction of the time series with MODIS images, but to a specific meteorological event. This anomaly was caused by meteorological anomalies in the short term. In particular, these disruptions observed in the RP of LST and JRP of LST and EVI/NDWI were mainly linked to the El Niño event that had a strong impact on the Amazonia in 2015 [35,51]. From the analysis of the RP of T max and H avg , other alterations can be noted from 2010 to 2014, which are not highlighted in the LST and even less so in the EVI and NDWI analysis. Since a single meteorological station was considered, this alteration was linked to local variations that cannot refer to the entire system.
To better identify the effect of drought events on the ecosystem function, we carried out the analysis using one-year sliding temporal windows on CRP, considering three years before and three years after 2005-2010 (AMO events) and 2015 (ENSO event). For each time window, we calculated the F-recurrence parameter for the couples EVI-NDWI, EVI-LST, and NDWI-LST. The F-recurrence varied between low and intermediate values, indicating low to intermediate synchronizations. Although the analysis explains why the values were relatively low, there were no statistically significant differences in the distribution of the F-recurrence between the period before and after the ENSO event. The p-values of the Student's T-test were always 5% higher than the threshold (Table 1) and, thus, we cannot reject the null hypothesis that the mean values of the F-recurrence distributions before and after the drought event are similar.
This suggests that the droughts did not produce a change in the synchronization of the pattern of the variables investigated (Tables 1-3). The F analysis for the EVI and NDWI indices was also performed for the four single sampling areas (located in different spatial locations) to verify the presence of a spatial different pattern caused by drought disturbance events. This analysis suggests that ENSO droughts events produced a change in the synchronization of the pattern of the investigated variables only for the area and the AMO events did not produce any effect (Table 4). The F analysis for the EVI and NDWI indices was also performed for the four single sampling 358 areas (located in different spatial locations) to verify the presence of a spatial different pattern caused 359 by drought disturbance events. This analysis suggests that ENSO droughts events produced a change 360 in the synchronization of the pattern of the investigated variables only for the area and the AMO 361 events did not produce any effect (Table 4).

372
The results obtained in this paper show that the methodology employed based on RQA is 373 capable of assessing the resilience of a natural system to a specific disturbance event. In fact, it 374 highlighted how the forest responded to a disturbing action due to drought events (global climate 375 effect), which did not alter the functionality of the forest.

Discussion
Several studies have reported on the impacts of earlier extreme temperatures and drought events on rainforests. With decreasing precipitation caused by drought events, the persistence of the structure and functions of the forest can decrease and shift to an alternative low tree-cover or low capacity to be a carbon sink [68]. However, the impacts of climate anomalies on ecosystem functions are still uncertain. Despite progress in recent years, the complex and non-linear vegetation-rainfall interactions are still poorly represented in process-based vegetation-climate studies [23,67,69].
The results obtained in this paper show that the methodology employed based on RQA is capable of assessing the resilience of a natural system to a specific disturbance event. In fact, it highlighted how the forest responded to a disturbing action due to drought events (global climate effect), which did not alter the functionality of the forest.
Specifically, the first part of the analysis showed the capacity of the area of the Amazon forest analyzed to maintain the same vegetation structure since the vegetation cover did not change at the large scale. At a smaller scale, this does not mean that drought events did not lead to the death of individual trees, because this is the main consequence of severe drought events [20,23,[25][26][27].
The analysis of the synchronization of several vegetation indices that are connected to photosynthesis represents a signature that characterizes a specific ecosystem. Variations in this signature imply variations in the system. Analyzing the functions of the forest with a recurrence analysis and considering the forest as a whole, the synchronization of the vegetation indices has shown that they were not altered by drought events. An alternative explanation is that the intensity of the extreme climate events was not sufficiently strong to produce changes in the recurrence pattern of the photosynthetically active green leaf canopy and water content.
Our study therefore suggests the good adaptation of the forest to seasonal drought at the ecosystem scale. This may be linked to the capacity of the roots to use the deep soil water, maintaining the water tension in the xylem relatively constant in both wet or dry seasons [70][71][72]. A combination of access to deep soil water, which helps to ensure evergreen canopies during dry seasons, and the low cloudiness during the dry season, helps the photosynthesis to continue [24]. This may also be linked to the higher quantity of light available for the forest during dry events, which is more limiting than water for tropical forest productivity [73,74].
Warm SSTs over the tropical Pacific (ENSO) induced changes in the atmospheric circulation, which had short temporal effects on the synchronization of the vegetation indices and are different from the changes induced by warm SSTs over the tropical Atlantic. The only variation detected was relative to the LST during ENSO events. This is influenced by the atmospheric temperature, but also by the capacity of the vegetation to mitigate the microclimate with the evapotranspiration of the plants. This local temporal alteration may thus also be interpreted as a lower ability of the forest to mitigate the local microclimate due to stress. However, the disruption was not irreversible and was highlighted only during the ENSO event. Thus, it did not affect the synchronization of the vegetation indices before and after the event.
Considering the similar structure of the vegetation forest, this suggests that the functionality of the forest in relation to photosynthetic activity may be similar before and after the event. The high capacity of the forest to absorb the disturbance can be explained by the feedback of the drought events which, aside from leading to the deaths of individual trees, can produce a "green up" at the ecosystem level, and which increases in solar radiation during drier periods that boost tropical productivity. Therefore, this creates a new opportunity to develop an ecological function with a higher resilience [75][76][77][78][79].
However, the results obtained in this study may be affected by the selection of sampling areas. Furthermore, the 1 km spatial resolution of the images and, especially, the spatially-averaged values, which might hide local variations, can reduce the capacity of the vegetation indices to discriminate the differences. Therefore, the results can be considered valid at the ecosystem scale. Moreover, the time series of spatial averages can hide the spatial variability that the drought events can generate. Repeating the analysis F between the EVI and NDWI for single samples at four different points revealed that only sample C had a different synchronization after the ENSO. This study can thus be used to identify both a global response of the forest system and specific individual local responses.
Focusing on the potential of the RQA, this study has however demonstrated the ability of the method to highlight different recurrence patterns of the vegetation time series caused by disturbance events.
Of course, this methodology does not replace other methodologies commonly employed in the literature, however, it should be considered as an approach that offers a different view in the study of the dynamics of natural systems.

Conclusions
Drought events generated by AMO in 2005 and 2010 did not induce functional changes in the forest ecosystem, while the drought event occurring in 2015, linked to ENSO, produced an immediate variation in LST, which, however, did not do any irreversible damage to the functionality of the forest system. The forest thus showed a good functional resilience to the drought events that occurred in the three years.
These findings refer to the specific RQA analysis applied to the selected ecosystem proxies. Specifically, the "retrospective analysis" using RQA applied to the vegetation indices enabled us to understand the past trajectory of a portion of the Amazon forest. The RQA highlighted the recurrent pattern of the natural processes linked to disturbance events in the short and medium term as well as the system response and its ability to recover.
In the first part of the analysis, the RQA evaluated the effect of the climatic disturbances on the evolution of the vegetation indices, that is, at the ecosystem level, the forest was able to overcome the extreme event of ENSO observed in 2015. Under continued global warming and a projected increase in the frequency of ENSO events, more frequent record-breaking climate extremes are expected in Amazonia [80]. We are thus confident that the approach followed here will be useful to evaluate the direct effects of these extreme events on the evolution of the capacity of the forest to maintain its ecological functions.
This approach could also be useful to understand the evolution of the system and its response after the fire that affected a large area of the forest in 2019. A change in the forest functions due to selective logging or below canopy burning, which are not easily detectable from greenness changes, could also be revealed using the RQA. This is important, because although of anthropogenic origin, fires in forested areas have been suggested as having the biggest impact on forest transition under an increased drought frequency scenario [63,81].
Overall, this work has highlighted the following potential of the RQA: • it provides an immediate and simple interpretation of the results through plot visualization; • it can be applied to short time series and allows an objective interpretation of specific variations; • it evaluates the ecological resilience, which relates to the magnitude of the disturbance that a system can absorb before shifting to another ecological condition [81]; and • it highlights the panarchy of processes at different spatial and temporal scales (e.g., the effect of the large scale ENSO on the smaller scale Amazon system).
Future developments could involve, for example, (i) extending the research to the entire Amazon forest to understand if the findings are representative of the whole biome; (ii) calibrations through field measurements that would improve the analysis; and (iii) the direct comparison of CO 2 emissions to better understand the feedback between human activities, climate conditions, and photosynthesis activity.
The vegetation indices approach clearly provides general information on the dynamics linked to the functions of the forest at the large scale, however, it cannot provide specific indications in relation to the small changes to individual plants caused by extreme events. Our analysis should therefore be integrated with local scale samplings.