Remote Sensing Phenology of the Brazilian Caatinga and Its Environmental Drivers

The Caatinga is the largest nucleus of Seasonally Dry Tropical Forests (SDTF) in the Neotropics. The leafing patterns of SDTF vegetation are adapted to the current environmental and climate variability, but the impacts of climate change tend to alter plants’ phenology. Thus, it is necessary to characterise phenological parameters and evaluate the relationship between vegetation and environmental drivers. From this information, it is possible to identify the dominant forces in the environment that trigger the phenological dynamics of the Caatinga. In this way, remote sensing represents an essential tool to investigate the phenology of vegetation, particularly as it has a long series of vegetation monitoring and allows relationships with different environmental drivers. This study has two objectives: (i) estimate phenological parameters using an Enhanced Vegetation Index (EVI) time-series over 20 years, and (ii) characterise the relationship between phenologic dynamics and environmental drivers. TIMESAT software was used to determine four phenological parameters: Start Of Season (SOS), End Of Season (EOS), Length Of Season (LOS), and Amplitude (AMPL). Boxplots, Pearson’s, and partial correlation coefficients defined relationships between phenologic dynamics and environmental drivers. The non-parametric test of Fligner–Killeen was used to test the interannual variability in SOS and EOS. Our results show that the seasonality of vegetation growth in the Caatinga was different in the three experimental sites. The SOS was the parameter that presented the greatest variability in the days of the year (DOY), reaching a variation of 117 days. The sites with the highest SOS variability are the same ones that showed the lowest EOS variation. In addition, the values of LOS and AMPL are directly linked to the annual distribution of rainfall, and the longer the rainy season, the greater their values are. The variability of the natural cycles of the environmental drivers that regulate the ecosystem’s phenology and the influence on the Caatinga’s natural dynamics indicated a greater sensitivity of the phenologic dynamics to water availability, with precipitation being the limiting factor of the phenologic dynamics. Highlights: The EVI time series was efficient in estimating phenological parameters. The high variability of the start of season (SOS) occurred in sites with low variability of end of the season (EOS) and vice versa. The precipitation and water deficit presented a higher correlation coefficient with phenological dynamics. Length of Season (LOS) and amplitude (AMPL) are directly linked to the annual distribution of rainfall.


Introduction
The Seasonally Dry Tropical Forests (SDTF) are characterised by a strong seasonal and irregular distribution of rainfall throughout the year, resulting in long dry season periods [1,2]. The Brazilian Caatinga, with an area of approximately 900,000 km 2 , is the largest nucleus of SDTF in the Neotropics [The tropical New World biogeographic region comprises Central America, the Caribbean, and South America] [3]. Moreover, with 3347 plant species, of which 526 are endemic, the Caatinga is the richest nucleus of SDTF in the Neotropics [4]. Its remarkable floristic diversity makes the Caatinga twice as rich as the Amazon rainforest when considering the species/area relationship [4]. Phenology studies recurring life-cycle events such as bird migration, flower blooming, or leaf emergence and senescence and the causes of their timing by biotic and abiotic forces [5,6]. The leafing patterns of Caatinga vegetation are adapted to the intense climate and water seasonality, being highly dependent on the climate interannual variability [7,8]. Most leaves fall during the dry season, and the first rainfall pulses trigger a quick leaf flush in the wet season [7][8][9].
Precipitation and soil moisture are the leading environmental drivers for the leaf changes in the Caatinga [8,9]. Still, the temperature can also be a driver for species presenting scheduled phenology in dry ecosystems [10]. The Intergovernmental Panel on Climate Change [11] forecasts an increase of 1.5 • C in the global mean air temperature for the next two decades in an optimistic scenario, considering the reduction of current emissions of CO 2 . Besides the temperature increases, climate changes are likely to alter the precipitation regimes in the following decades [12,13]. The consequences of these changes are, among others, the projected drying out of surface soils [12] and prolonged dry seasons, with an increase of 47% of the area of the Northeast Region of Brazil (NEB) subjected to extreme drought events until 2070 [13]. Thus, understanding the past and current vegetation's response to the environmental drivers is paramount to predicting its behaviour in climate change scenarios, allowing the detection of changes in the timing of leaf patterns and their causes.
The ability to monitor global vegetation phenology, or Land Surface Phenology (LSP), has increased with the validated Remote Sensing (RS) and modelling approaches to mapping phenology [14,15]. Long-term data from satellite products are useful tools for understanding the phenological responses of vegetation to current environmental drivers using Vegetation Indices (VIs), allowing it to predict its responses to climate change scenarios. VIs time series has received the attention given its potential to characterise interactions between climate and vegetation with broad applications in different ecosystems [16][17][18]. Several VIs are calculated based upon different spectral bands and, therefore, evidence of different components of the environment [19]. The Enhanced Vegetation Index (EVI) has been widely used to characterise vegetation phenology [20,21] due to its sensitivity to high biomass and reduced atmospheric and soil effects. EVI is calculated from the near-infrared (NIR), Red, and Blue bands and can be derived for different satellite platforms, such as Landsat, Sentinel, and MODIS. The use of algorithms to determine the main phenological metrics from the VIs time series has favoured the representation of the phenological stages of each cropping system, allowing a crop-type classification based on their phenological metrics [22]. However, the studies driving this on a global scale have been primarily focused on forest ecosystems, associating phenological changes in vegetation with climate patterns, particularly with rainfall data [23]. These studies showed that the phenological dynamics strongly depend on the seasonality of rainfall [21,24]. Still, the studies on a regional scale indicate that other environmental drivers also trigger phenological changes [25][26][27][28].
The LSP applied to an ecosystem scale seems to offer the best opportunities to advance understanding of environmental triggers and determinants for phenological dynamics, given the possibility to understand it on a broad scale, encompassing areas in a range of contrasting environmental conditions. For instance, the early greening or pre-rain green-up, a phenomenon where trees produce leaves before the rain starts, was registered in the woodlands and savannas of southern Africa through RS satellite techniques [29,30]. Furthermore, the application of the LSP at a continental scale and using long-term time-series  (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) allowed us to measure the variability in leaf flushing (i.e., greening) among years and to identify the photoperiod as the environmental cue for early greening [30].
The use of LSP and their drivers will be significant for SDTFs where interannual rainfall variability and rainy season duration change on a spatial and temporal scale [31,32], factors that are expected to influence the phenological strategies of plant communities in this vegetation [10]. There was also great regional variability and interannual fluctuation in vegetation phenology, and the overall phenological trends shifted later [21]. By following the studies for dry forests, Tong et al. [15] reported that the interannual rainfall variability was a more dominant force than fire events and land-use change in the phenological trend in tropical areas. For example, Jesus et al. [25] noticed changes in the phenological patterns for dense and open vegetation areas of the Caatinga, suggesting that factors that vary at spatial scales, such as the vegetation structure, would also be necessary for the phenological responses of the vegetation. In an experimental area in the Caatinga, the phenological response was directly related to soil water availability [9].
Despite advances in the analysis of phenological patterns and their associated environmental drivers, mainly observed in studies at the ecosystem scale, the application of long-term time series of vegetation indices in studies of the Caatinga vegetation is scarce. In addition, there are limitations to the diversity of sites studied for the Caatinga. For example, when analysing the dominant environmental drivers for the phenology of seasonally dry ecosystems (Caatinga, Cerrado), Alberton et al. [33] observed that the dominant drivers in these ecosystems were distinct, with light (measured as day-length) being more relevant in explaining leafing patterns in Cerrado communities than rainfall for Caatinga communities. Therefore, comparing sites of the same ecosystem can better define the environmental drivers associated with the phenological dynamics. There are also limitations to the number of environmental drivers analysed in the studies carried out for the Caatinga. Analysis with more environmental drivers could reveal meaningful soil-plant-atmosphere interactions, which may occur to a lesser extent.
Given the above, this study proposes to evaluate the Caatinga phenological sensitivity to environmental drivers in three Caatinga ecoregions. The seasonality of vegetation will be observed from the EVI time series over 20 years (2000-2019) and environmental drivers (precipitation, air temperature, soil moisture, and water deficit) from global databases. This study has two objectives: (i) estimate phenological parameters using an EVI time-series over 20 years, and (ii) characterise the relationship between phenologic dynamics and environmental drivers. The results will also be expected to serve as a baseline against which to compare future changes in Caatinga phenology due to natural or anthropogenic causes.

Study Areas
The semi-arid climate (Köppen's BSh) [34], where the Caatinga is located, presents accumulated precipitation below 600 mm and an average air temperature of 25 • C to 30 • C throughout the year [35]. However, seasonally, the average annual rainfall varies widely from 300 mm to 1000 mm annually, mainly concentrated in a period of 3 to 4 months, during summer and autumn, followed by a prolonged dry season lasting 8 to 9 months during winter and spring [36]. The Caatinga is characterised by high interannual precipitation variability, with droughts that can last for years and high potential evapotranspiration rates of between 1500 mm and 2000 mm annually [37]. The Caatinga region has different soil categories, ranging from shallow, rocky, and relatively fertile to deep with high natural fertility and sedimentary or sandy, deep, and non-fertile [35]. In the Caatinga, the typical vegetation comprises deciduous and thorny forests or small forests that grow mainly in exposed crystalline rocky terrains, such as the Depression Sertaneja [38]. The region's geomorphology ranges from lowlands of 300 m, mountains, plateaus, and plateaus with an altitude of up to 1000 m [35].
The study was carried out at three SDTF protected sites, each area being in a different ecoregion (   [39]. The average annual temperature is 26.2 • C, and the rainfall is 387 mm [39], concentrated mainly from January to April. The predominant soil type in the experimental site is Red-Yellow Argisol, followed by Haplossol and Hydromorphic Vertisol [40]. The SD experimental site's geographic coordinates are 6 • 34 42 S, 37 • 15 05 W, and 203 m altitude. Its area is characterised by dry xerophilous forest with deciduous plant species and the predominance of small, widely dispersed trees and shrubs up to 10 metres high, which develop and grow only in the rainy season between January and May [41]. The average annual precipitation varies between 400 and 700 mm, and the average yearly temperature ranges from 28 to 30 • C. The dominant soil is Neosol Litolic, shallow (about 40 cm), stony, and low fertility [42]. The BP experimental site's geographic coordinates are 7 • 16 49 S, 35 • 58 34 W, and 492 m altitude. The vegetation present in the experimental site is composed of dense Caatinga, with tree heights between three and seven metres. The region experiences a rainy season starting in March/April and extending until August, while the dry season lasts for five to six months [43]. The average annual rainfall of Campina Grande is about 750 mm, and the average annual air temperature is 23.3 • C [44]. The soil is nitric planosol [45].
forests that grow mainly in exposed crystalline rocky terrains, such as the Depression Sertaneja [38]. The region's geomorphology ranges from lowlands of 300 m, mountains, plateaus, and plateaus with an altitude of up to 1000 m [35].
The study was carried out at three SDTF protected sites, each area being in a different ecoregion ( Figure 1): (i) Meridional Sertaneja Depression (MSD) ecoregion experimental site located Embrapa Semi-arid, Petrolina, Pernambuco; (ii) Sertanejas Depressions (SD) ecoregion experimental site located at Seridó Ecological Station, Serra Negra do Norte, Paraíba; (iii) Borborema Plateau (BP) ecoregion experimental site located at Semi-arid National Institute, Campina Grande, Paraíba. The MSD experimental site's geographic coordinates are 9°02′47″S, 40°19′16″W, 395 m altitude, and its vegetation is dry xerophilous forest, with trees and shrubs that average five metres in height [39]. The average annual temperature is 26.2 °C, and the rainfall is 387 mm [39], concentrated mainly from January to April. The predominant soil type in the experimental site is Red-Yellow Argisol, followed by Haplossol and Hydromorphic Vertisol [40]. The SD experimental site's geographic coordinates are 6°34′42″S, 37°15′05″W, and 203 m altitude. Its area is characterised by dry xerophilous forest with deciduous plant species and the predominance of small, widely dispersed trees and shrubs up to 10 metres high, which develop and grow only in the rainy season between January and May [41]. The average annual precipitation varies between 400 and 700 mm, and the average yearly temperature ranges from 28 to 30 °C. The dominant soil is Neosol Litolic, shallow (about 40 cm), stony, and low fertility [42]. The BP experimental site's geographic coordinates are 7°16′49″S, 35°58′34″W, and 492 m altitude. The vegetation present in the experimental site is composed of dense Caatinga, with tree heights between three and seven metres. The region experiences a rainy season starting in March/April and extending until August, while the dry season lasts for five to six months [43]. The average annual rainfall of Campina Grande is about 750 mm, and the average annual air temperature is 23.3 °C [44]. The soil is nitric planosol [45].

Data Processing
The processes for obtaining geospatial images and data were performed using the Google Earth Engine (GEE) tool [46]. Cloud computing platforms, such as GEE, facilitate the processing of satellite images globally. With access and analysis of several petabytes of images and a broad set of geospatial data, dating from 40 years of Earth observation data, updated and expanded daily, the tool allows viewing, manipulating, editing, and creating spatial data [46]. The EVI vegetation indices derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor reflectance data were used to represent the variability of native vegetation, investigating some biophysical mechanisms that involve soil-plant-atmosphere relationships in the natural dynamics of the Caatinga. The MODIS sensor images and the TerraClimate and CHIRPS data sets' environmental drivers were evaluated from 2000 to 2019. The environmental drivers were made up of monthly precipitation, air temperature, soil moisture, and water deficit time series for the three experimental areas.

MODIS Data
This research used the surface reflectance data of the MODIS 16-day nadir BRDF-Adjusted Reflectance product (MCD43A4); the dataset is produced daily using 16 days of Terra and Aqua MODIS data at 500 metres resolution. The reflectance bands were used to obtain the EVI vegetation index [47]. The first available images of MODIS date from February 2000, and this study had its period limited to including images until December 2019. The EVI series were assembled for each studied site, with the value of the MCD43A4 product pixel corresponding to the geographic coordinates of each testing location. After that, the EVI time series was converted into a monthly EVI time series, using the mean value of each month.

TerraClimate e CHIRPS Data
The air temperature, precipitation, soil moisture, and water deficit are some of the potential environmental drivers for the leaf exchange strategies of Caatinga plant species [33]. These possible environmental drivers acting on the Caatinga vegetation were observed in three different sites of the Caatinga. All the variables mentioned were obtained monthly from the TerraClimate data set [48], except for precipitation data. Monthly precipitation data were obtained from the Climate Hazards Group InfraRed Precipitation with Stations (CHIRPS) database [49]. The exact centre pixels used to construct the EVI time series were used to obtain the environmental drivers' time series.

Phenological Metrics
The phenological metrics represent the characteristics of the vegetation within its phenological cycle, or phenophases, corresponding to dimensionless output parameters and can be calculated based on the EVI time series. In this study, the TIMESAT software [50,51] was used to analyse 20 years of EVI time series (from 2000 to 2019) and to compute 4 phenological metrics: Start of Season (SOS), End Of Season (EOS), Length of the season (LOS), and Amplitude (AMPL) difference between the peak and the base level value. After applying the Savitzky-Golay filter [52] in TIMESAT, a seasonality parameter per year was chosen, representing a phenological cycle with a start and end level of 20% of the seasonal amplitude. This threshold value was used in several studies and is known to be accurate in registering the plant's phenological transitions [53][54][55][56][57][58]. The four phenological metrics are shown in boxplot graphics for each ecoregion studied and used for correlations with environmental drivers.

Seasonal Variability Analysis
Interannual variability in phenological metrics across sites-To test if the interannual variability in SOS and EOS (i.e., the variances in phenological transition dates) differs among the three sites, we used the non-parametric test Fligner-Killeen [59]. The Fligner-Killeen test compares the homogeneity of variances among samples [59]. The same analysis was performed to test if the variability of SOS differs from the EOS within each site. To test if the LOS and the AMPL of EVI differ among sites, we performed a one-way analysis of variance (ANOVA) followed by the Tukey post-hoc test. In addition, we used box plots to compare the three experimental sites' seasonal data on EVI, phenological metrics, and environmental drivers for 20 years.
Influence of environmental drivers on phenological metrics-The environmental drivers were also observed monthly to assess their influence on phenological parameters. This study analysed these relationships through scatter plots with smoothed (r) correlation curves from Pearson's classification (P) between the EVI and environmental drivers. Pearson's correlation between EVI time series and monthly environmental drivers was submitted to different monthly lag periods (0-3). The lag monthly of each environmental driver with the highest Pearson's correlation was applied to the partial correlation method. Then the partial correlation method was used to analyse the environmental drivers in the monthly EVI time series [60]. Data analysis was carried out in R [61].
The monthly values of the seasonal EVI profile for each ecoregion studied using 20 years of time series are presented in the box plot of Figure 2 Although the graphs with the monthly EVI values presented in Figure 2 understand the temporal behaviour, they do not reveal as many characteristics about the vegetation as the graphs with the phenological metrics presented in Figure 3. When analysing the SOS, it is noticed that the BP and MSD experimental sites present a greater interquartile range when compared to the SD experimental site. For most observations, the SOS from the SD experimental site has taken place between January and February. On the MSD experimental site, it is observed that SOS occurs most frequently between November and January. At BP, it is observed that the highest frequency of SOS occurs between January and May (117 days).
Interestingly, the interquartile range observed for the SOS is not for the EOS and LOS metrics for the MSD and BP experimental sites. For the EOS and LOS, it was observed that the MSD and BP sites have a smaller interquartile range than the SD experimental site, with interquartile range values of less than 30 days for both situations. For the SD experimental site, the variation in the interquartile range in the EOS is similar to that observed for SOS, with a variation slightly greater than 30 days interval and the LOS reaching an interquartile range of 60 days. A relationship of the amplitude with the other metrics, SOS, EOS, and LOS, was not observed. The amplitude's highest values were observed for the BP and SD experimental sites, with 25% of the observations above 0. 46 Although the graphs with the monthly EVI values presented in Figure 2 understand the temporal behaviour, they do not reveal as many characteristics about the vegetation as the graphs with the phenological metrics presented in Figure 3. When analysing the SOS, it is noticed that the BP and MSD experimental sites present a greater interquartile range when compared to the SD experimental site. For most observations, the SOS from the SD experimental site has taken place between January and February. On the MSD experimental site, it is observed that SOS occurs most frequently between November and January. At BP, it is observed that the highest frequency of SOS occurs between January and May (117 days).  Interestingly, the interquartile range observed for the SOS is not for the EOS and LOS metrics for the MSD and BP experimental sites. For the EOS and LOS, it was observed tha the MSD and BP sites have a smaller interquartile range than the SD experimental site with interquartile range values of less than 30 days for both situations. For the SD

Environmental Drivers
Due to the high number of deciduous species, the SDTF, such as the Caatinga, presents high variability in plant biomass in the annual cycle. Figures 4-6 present the environmental drivers' boxplots that may influence the triggers of phenological changes in the Caatinga. When observing the shape of the graphs of the median behaviour of each environmental driver, it can be seen that there is a similarity in the water deficit, precipitation, and soil moisture behaviour. Pearson's correlation coefficient for the environmental drivers and time series EVI is presented in Figure 7. Figure 7 only presents the lag (0-3 months) environmental drivers with the highest Pearson's correlation. With the application of Pearson's correlation, this study allowed for analysis of each environmental driver and identified their action times on vegetation. These results reinforce the importance of precipitation as the environmental driver that best reflects plant biomass production in the Caatinga. Thus, it is the environmental driver that presents the highest coefficient correlation among all: MSD (r = 0.7258; p < 0.05, lag = 1), SD (r = 0.8267; p < 0.05, lag = 1), and BP (r = 0.7546; p < 0.05, lag = 1). Subsequently, water deficit and soil moisture had the highest correlation values. The water deficit's correlations were: MSD (r = −0.6, p < 0.01, lag = 0); SD (r = −0.79, p < 0.01, lag = 0); BP (r = −0.69, p < 0.01, lag = 1). Soil moisture is one of the main environmental triggers of the Caatinga, and in shallow soils, there is a tendency to saturate and also dry out more quickly, not allowing, in many cases, a long-term response from the vegetation: MSD (r = 0.52; p < 0.05, lag = 0), SD (r = 0.69; p < 0.05, lag = 0) and BP (r = 0.52; p < 0.05, lag = 0). Pearson's correlation coefficient between the EVI and air temperature time series showed the weakest but most significant correlation. While MSD and SD presented positive relations between EVI and temperature (r = 0.6; p < 0.05, lag = 2) and (r = 0.59; p < 0.05, lag = 3), BP (r = −0.56; p < 0.05, lag = 0) showed a negative association.         In Pearson's correlation, environmental drivers were analysed one by one. However, the observed behaviour may not reflect the actual effects on vegetation represented by the EVI time series. There may be an association of impact between the analysed environmental drivers. Partial correlation analysis is a way to solve this problem, allowing the analysis of multiple variables. The partial correlation analysis between the environmental drivers and EVI is presented in Table 1 with five scenarios. Scenario 1 considers all environmental drivers to calculate partial correlations, and in the other scenarios, the effect of one of the drivers is retained. Precipitation and soil moisture always showed a positive partial correlation. In contrast, the water deficit showed negative correlations, with r values ranging from −0.25 to −0.62. Only at the BP site, the partial correlation with temperature was negative. In scenario 1, with all environmental drivers, precipitation was the environmental driver that presented the highest correlation, with r values ranging from 0.48 to 0.51. The other environmental drivers are better perceived In Pearson's correlation, environmental drivers were analysed one by one. However, the observed behaviour may not reflect the actual effects on vegetation represented by the EVI time series. There may be an association of impact between the analysed environmental drivers. Partial correlation analysis is a way to solve this problem, allowing the analysis of multiple variables. The partial correlation analysis between the environmental drivers and EVI is presented in Table 1 with five scenarios. Scenario 1 considers all environmental drivers to calculate partial correlations, and in the other scenarios, the effect of one of the drivers is retained. Precipitation and soil moisture always showed a positive partial correlation. In contrast, the water deficit showed negative correlations, with r values ranging from −0.25 to −0.62. Only at the BP site, the partial correlation with temperature was negative. In scenario 1, with all environmental drivers, precipitation was the environmental driver that presented the highest correlation, with r values ranging from 0.48 to 0.51. The other environmental drivers are better perceived when precipitation is removed from the analysis (scenario 2). In scenario 2, the highest correlations were for temperature in MSD (0.46), water deficit and soil moisture for SD (0.40), and water deficit for BP (−0.62). This dependence relationship becomes evident when the drivers with the highest partial correlation identified in scenario 2 are removed from the analysis. Thus, the highest values of r for precipitation are seen in scenarios 5, 3, and 4 for MSD, SD, and BP sites, respectively. Table 1. Partial correlation coefficients between the EVI and environment drivers across MSD, SD, and BP sites. Scenario 1 considers all environmental drivers to calculate partial correlations. For the other scenarios, the effect of one of the variables is removed: scenarios 2 (without precipitation), 3 (without soil moisture), 4 (without water deficit), and 5 (without temperature). "-" indicates that the environment driver was not used to calculate the partial correlation. Only for statistical significance (p < 0.05). Missing estimates (NS) are not significant.

Discussion
Compared to the following season, the greater rainfall variability at the beginning of the wet season may be the primary driver for the interquartile range of the SOS to be greater than that observed for the EOS and LOS (Figure 3). The higher frequency and accumulation are the main characteristics of rainfall contributing to the higher AMPL and LOS values at the BP experimental site. The partial coefficient of air temperature in the BP experimental site was always lower. In tropical regions, the air temperature variability is low, with about four degrees of annual amplitude difference from the monthly normal. However, in the MSD experimental site, there was a greater partial correlation between air temperature (scenario 2). The climatological water deficit is calculated from the reference evapotranspiration and precipitation, providing an efficient measure of the water availability and demand of the environment. Air temperature, directly related to water vapour-pressure deficit, is one of the climatic factors contributing to the variability of reference evapotranspiration [62]. The average monthly precipitation is quite variable in the annual cycle, which has a greater impact on the water deficit.
In the Caatinga, vegetation's phenology and photosynthetic activity are associated with water availability [9]. In drier regions such as MSD and SD, intra-annual water availability is lower, leading to shorter growing season periods. The water stress makes the environmental conditions (Figures 2 and 3) harsh for leaf maintenance for vegetation. In contrast, BP with more regular water availability (Figures 2 and 3) shows larger LOS, which means more water and carbon exchange time at the soil-vegetation-atmosphere interface [63]. Carbon assimilation was related to EVI at Caatinga by Mendes et al. [64], demonstrating that the vegetation presents more significant photosynthetic activity and productivity (seasonal and total) during high EVI value periods. The vegetation at BP showed a higher senescence rate than in the other areas. The phenological parameter is related to leaf loss in the dry season and the maximum use of water availability for its recovery after the dry period with accelerated regrowth and increased metabolic processes in the rainy season [65]. These vegetation strategies are adaptive to optimise the phenological, vegetative, and reproductive processes [10]. The study areas showed similarities in the seasonality of phenological parameters, corroborating the high association level with water availability.
Many studies have shown that rainfall seasonality regulates SDTF canopy seasonality [24,66], but there is a complex relationship between environmental drivers and the vegetation response [9,28], mainly for the Caatinga where the plant physiology is adapted to drought and elevated temperatures [67]. As expected for the Caatinga vegetation, the months with the highest water deficit have the lowest leaf cover (lowest EVI values), a result previously found by other studies in the region [7,33,68]. About 70% of the year has a water deficit [69]. According to Flerchinger et al. [70], about 90% of the rainwater in arid and semi-arid regions returns to the atmosphere through evapotranspiration. As there is not such a significant variability in the average air temperature between the tropics, the annual constancy provides a high evaporative demand from the atmosphere throughout the year [69], conditioning the water deficit to the seasonal fluctuation of the rainfall, adjusting the Caatinga phenological cycle to the water availability. Phenological transitions are an excellent indicator of climate change [71], and future scenarios estimate greater water demand from plants and the occurrence of droughts, providing a more significant water deficit [72,73], and with that, being able to alter the phenology of the plants. The water deficit slows down plants' growth, causes leaves, fruits, and flowers to fall, and, in the short term, tends to anticipate flowering and the beginning of fruiting, reducing plant cycles. In contrast, they tend to extend or even prevent the regular completion of the plant cycle [74].
The months with the highest spectral response of vegetation occurred in the rainy season when rain and soil moisture were predominantly distributed. The peak of the EVI was preceded by the month of greatest precipitation at SD and MSD. The rain at the BP experimental site shows a smoother and more regular distribution than in other areas, with the EVI closely following seasonal fluctuations in precipitation. The ecosystem accumulates sufficient water reserves in the soil and biomass for both sites under study, resulting in slower leaf fall during the dry season. The EVI followed the monthly rainfall distributions linearly, as shown in Figure 4. Likewise, the air temperature is observed after the maximum peak of the EVI. The air temperature also decreases the vegetation response, resulting in greater water stress for the research areas. Because of the high diversity of species at the Caatinga, the studied sites could present plants at different stages of adaptation or with physiological aspects related to water uptake or leaf abscission. The fact that MSD and SD present shorter LOS than BP could be related to vegetation strategies to use the water more efficiently during the available period. The peak of EVI values was higher at MSD and SD, and this would be related to more biomass production resulting from the water use efficiency (WUE) of the species composition.

Conclusions
This study observed that using the EVI time series efficiently estimated phenological parameters. The determination of the phenological parameters for 20 years allowed the evaluation of their seasonality. The seasonality of Caatinga vegetation growth differed between the three experimental sites. The SOS was the parameter that presented the greatest variability in the days of the year (DOY), reaching a variation of 117 days. The high variability of SOS occurred in sites with low variability of EOS and vice versa. The relationship between air temperature and vegetation was insignificant for the two experimental sites and presented a significant value for the site drier (MSD). At the experimental site where the highest vegetative peak co-occurred with the low air temperature, the highest LOS and AMPL were observed. In addition, the values of LOS and AMPL are directly linked to the annual distribution of rainfall, and the longer the rainy season, the greater their values are. The variability of the natural cycles of the environmental drivers that regulate the ecosystem's phenology and the influence on the Caatinga's natural dynamics indicated a greater sensitivity of the vegetation to water availability, with the water deficit being the limiting factor of the vegetation. There is a greater need for new research and studies aiming to reduce the impact of climatic risks by characterising and quantifying species' response patterns to likely climate change. There will be changes in phenological patterns, production, and spatial distribution. Quantifying these impacts and devising appropriate actions to minimise such adversities represent significant research challenges in the coming decades.