Normalized Difference Vegetation Index Temporal Responses to Temperature and Precipitation in Arid Rangelands

Rangeland degradation caused by increasing misuses remains a global concern. Rangelands have a remarkable spatiotemporal heterogeneity, making them suitable to be monitored with remote sensing. Among the remotely sensed vegetation indices, Normalized Difference Vegetation Index (NDVI) is most used in ecology and agriculture. In this paper, we research the relationship of NDVI with temperature, precipitation, and Aridity Index (AI) in four different arid rangeland areas in Spain’s southeast. We focus on the interphase variability, studying time series from 2002 to 2019 with regression analysis and lagged correlation at two different spatial resolutions (500 × 500 and 250 × 250 m2) to understand NDVI response to meteorological variables. Intraseasonal phases were defined based on NDVI patterns. Strong correlation with temperature was reported in phases with high precipitations. The correlation between NDVI and meteorological series showed a time lag effect depending on the area, phase, and variable observed. Differences were found between the two resolutions, showing a stronger relationship with the finer one. Land uses and management affected the NDVI dynamics heavily strongly linked to temperature and water availability. The relationship between AI and NDVI clustered the areas in two groups. The intraphases variability is a crucial aspect of NDVI dynamics, particularly in arid regions.


Introduction
Rangelands cover almost 33% of ice-free land globally. Monitoring rangeland is a key aspect of stopping its degradation. Severe degradation is causing the disappearance of 5-10 million hectares of agricultural lands every year [1,2]. Rangelands management is complex as rangelands are not spatially homogeneous; those with erodible soils and palatable vegetation are more prone to degrade than others. Arid and semiarid rangelands tend to have erodible soils and variable precipitation regimes, making them more sensitive to degradation [3,4]. In addition, different types of rangeland depend on land use, plant communities, soil, and climatic conditions. These include grazed wastelands such as stubble from mainly rainfed cereal crops, other croplands and fallow lands, grasslands, scrublands, and open woodlands [5]. Using remote sensing to monitor these areas and improving the land's management practices can aid in preventing the degradation and reinforcing the sustainability of these ecosystems [2].
Remote sensing techniques are increasingly used by ecologists and agriculturalists [6][7][8][9], as a cost-effective means to measure the characteristics, biomass, and extent of vegetation [10]. Normalized Difference Vegetation Index (NDVI) is the most widely used vegetation index by rangeland ecologists [11]. It can be obtained at different spatial resolutions, the most common at 500 × 500 m 2 and 250 × 250 m 2 , named Low and Medium Resolution (LR and MR). NDVI shows a good correlation with biomass of different vegetation types in arid and semiarid areas [12,13]. Another useful index is the NDVI anomaly (Z NDVI ). This index has shown robust results to identify damaged vegetation, and therefore it is especially suited to analyze the status of semiarid and arid areas [14].
Temperature and precipitation are mostly studied as drivers of NDVI [15][16][17][18][19][20][21][22][23]. However, the interactions among them and with other factors are still not fully comprehended [24]. Temperature effects on NDVI are significant when water is available in the ecosystem; in these circumstances, precipitation plays a minor role. However, this relationship grows more assertive in arid regions, where water availability seems to be one of the main drivers of NDVI [25][26][27]. Furthermore, some authors [28] suggest studying additional climatic factors further to understand thermal and hydric stress effects on NDVI patterns. Other climatic variables that have been reviewed are soil moisture, evapotranspiration, and land use cover [16,22,23]. The relationships of NDVI with these climatic variables differ by season and vegetation type and biome. Another factor that makes it more difficult to disentangle the interactions is the effect of human activities on both NDVI and the ecosystem itself, especially when it comes to overexploitation such as overgrazing or water overuse [29]. Therefore, it is essential to deeply understand NDVI and meteorological variables' relationship, especially accounting for the differences in this relationship between seasons and across types of land management. This is particularly relevant for agrometeorological indexes that should be calculated during the year, as it is in the case of rangelands [30].
NDVI and its relationship with meteorological variables have reported different results depending on the analyzed spatial scale. Tarnavsky et al. (2008) [31] in Arizona (USA) and Stefanov et al. (2005) [32] in different locations worldwide, described differences in spatial heterogeneity. Peng et al. (2017) [33] studied different spatial scales in China, which exhibited different correlations depending on coarser or finer scales; all were performed at monthly or annual temporal scale. On the other hand, human activities negatively correlated with NDVI at a Medium Resolution (MR) while it exerted a positive correlation at a lower resolution. Therefore, considering different spatial scales to identify the causes of NDVI patterns remains a challenging task.
The Aridity Index (AI) represents water availability, and different expressions have been developed [34]. We used a modified version of a widely accepted ratio of annual precipitation to the annual potential evapotranspiration in this work, developed by the United Nations Environment Programme in 1992 [35,36]. This index has been used to quantify droughts and estimate possible changes in climate regimes [37,38], and manage afforestation and reforestation projects and prioritize and assess future conservation efforts in rangelands [39,40]. It has also been previously compared to vegetation indexes [41].
The objective of the study was to evaluate the responses of the rangeland NDVI to temporal dynamics of temperature and precipitation for an arid environment. To evaluate it, NDVI series from four areas, corresponding to different types of rangeland in Murcia (southeast of Spain), were acquired for 2002 to 2019, at two spatial resolutions (LR and MR). The correlation between the NDVI and the climatic variables and the AI was used to gain better insight into the patterns of rangelands to monitor and characterize them and optimize the management accordingly.

Area of Study
The area selected is located in Murcia province, in the southeast of Spain. The area has a Mediterranean arid climate with annual precipitation less than 300 mm, although there are regional variations [42]. Four square areas with 2.5 km sides were selected for this study. They are situated in the vicinity of a meteorological station, covering three different agricultural regions of Murcia: Northeast, Segura River, and Northwest ( Figure 1, Table A1 in Appendix A). The average temperature varies amongst the four areas from 14.7 to 17.3 • C, and the average accumulated precipitation per eight-day period ranges from 262 to 348 mm. emote Sens. 2021, 13, x FOR PEER REVIEW

Area of Study
The area selected is located in Murcia province, in the southeast of Sp has a Mediterranean arid climate with annual precipitation less than 300 m there are regional variations [42]. Four square areas with 2.5 km sides wer this study. They are situated in the vicinity of a meteorological station, c different agricultural regions of Murcia: Northeast, Segura River, and Nort 1, Table A1 in Appendix A). The average temperature varies amongst the fo 14.7 to 17.3 °C, and the average accumulated precipitation per eight-day p from 262 to 348 mm. All the areas are used as rangeland, two predominantly herbaceous (A1 two mainly covered by trees (A3 and A4). Area 1 (A1) is mostly covered by cereal crops. Area 2 (A2) is almost entirely covered by mixed croplands us grazing with some grassland and scrubland. Area 3 (A3) has a top grass with shrubs and with few forested areas surrounded by tree crops with irr cerning Area 4 (A4), it is mainly covered by coniferous woodland with m small patches ( Figure 2). For the following analyses, 11 and 61 pixels corr rainfed areas were selected for A3 LR and MR, respectively, to eliminate th rigation in the NDVI dynamics. The other areas are mainly conformed by grasslands, and reforestation that do not present irrigation [43]. Therefore, a were used since no irrigation could affect their NDVI dynamics. The selecti based on each pixel's average NDVI for summer months (June, July, and A with a summer average below 30 and did not have peaks over 40 in its ori ries, were selected to be analyzed, and these pixels match the grassy patch th from the top center to the right bottom corner. All the areas are used as rangeland, two predominantly herbaceous (A1 and A2) and two mainly covered by trees (A3 and A4). Area 1 (A1) is mostly covered by stubble from cereal crops. Area 2 (A2) is almost entirely covered by mixed croplands used for stubble grazing with some grassland and scrubland. Area 3 (A3) has a top grassy area mixed with shrubs and with few forested areas surrounded by tree crops with irrigation. Concerning Area 4 (A4), it is mainly covered by coniferous woodland with mixed crops on small patches ( Figure 2). For the following analyses, 11 and 61 pixels corresponding to rainfed areas were selected for A3 LR and MR, respectively, to eliminate the effects of irrigation in the NDVI dynamics. The other areas are mainly conformed by rainfed crops, grasslands, and reforestation that do not present irrigation [43]. Therefore, all their pixels were used since no irrigation could affect their NDVI dynamics. The selection was made based on each pixel's average NDVI for summer months (June, July, and August). Pixels with a summer average below 30 and did not have peaks over 40 in its original time series, were selected to be analyzed, and these pixels match the grassy patch that crosses A3 from the top center to the right bottom corner.

NDVI and Meteorological Data Collection
To study spatial resolution effects, the target areas were studied using 500 × 500 m 2 , 25 pixels for each area, and 250 × 250 m 2 , implying a set of 132 pixels ( Figure 2). MOD09Q1.006 (LR) and MOD09A1.006 (MR) were collected from AppEEARS [44], downloading the RED (band 1) and NIR (band 2) values for the target areas at these spatial scales. Both had an 8-day temporal resolution from the beginning of 2002 through 2019, a total of 18 years of data. A maximum temporal resolution among the cited literature was found to be biweekly [16], and more commonly by month [15,[17][18][19][20][21]. For each pixel series, R was used to calculate the Normalized Difference Vegetation Index (NDVI), using the following formula: The possible NDVI values range from 0 to 100. The NDVI obtained values were then checked for quality. If the data were not categorized as ideal quality in the quality band from AppEEARS, less than 0.01% was deleted for every area. The gaps were filled using running averages with a gap interval of seven dates. The time series were then smoothed using the Savitzky-Golay method [45], with a window size of nine selected, based on the best-fitted outputs. The NDVI anomaly (ZNDVI) was calculated by applying a z-score by date [46]: where µNDVI is the average of all data in all available years that were measured in the same calendar date and σ is the corresponding standard deviation. In this new series, the seasonal variations of NDVI were eliminated. Daily meteorological data from the closest meteorological stations (Table A1) were also used [47]. Average temperature and accumulated precipitation were calculated for every eight days to match the NDVI dates.
To examine the variation of NDVI during the year, boxplots of NDVI and meteorological variables were used. Different phases were defined based on NDVI pattern during the year as this region presents harshness and aridity. These phases ( ) were based on the trend of NDVI values: increasing, decreasing, or constant. A Chow test [48] of NDVI

NDVI and Meteorological Data Collection
To study spatial resolution effects, the target areas were studied using 500 × 500 m 2 , 25 pixels for each area, and 250 × 250 m 2 , implying a set of 132 pixels ( Figure 2). MOD09Q1.006 (LR) and MOD09A1.006 (MR) were collected from AppEEARS [44], downloading the RED (band 1) and NIR (band 2) values for the target areas at these spatial scales. Both had an 8-day temporal resolution from the beginning of 2002 through 2019, a total of 18 years of data. A maximum temporal resolution among the cited literature was found to be biweekly [16], and more commonly by month [15,[17][18][19][20][21]. For each pixel series, R was used to calculate the Normalized Difference Vegetation Index (NDVI), using the following formula: The possible NDVI values range from 0 to 100. The NDVI obtained values were then checked for quality. If the data were not categorized as ideal quality in the quality band from AppEEARS, less than 0.01% was deleted for every area. The gaps were filled using running averages with a gap interval of seven dates. The time series were then smoothed using the Savitzky-Golay method [45], with a window size of nine selected, based on the best-fitted outputs. The NDVI anomaly (Z NDVI ) was calculated by applying a z-score by date [46]: where µ NDVI is the average of all data in all available years that were measured in the same calendar date and σ NDVI is the corresponding standard deviation. In this new series, the seasonal variations of NDVI were eliminated. Daily meteorological data from the closest meteorological stations (Table A1) were also used [47]. Average temperature and accumulated precipitation were calculated for every eight days to match the NDVI dates.
To examine the variation of NDVI during the year, boxplots of NDVI and meteorological variables were used. Different phases were defined based on NDVI pattern during the year as this region presents harshness and aridity. These phases (i) were based on the trend of NDVI values: increasing, decreasing, or constant. A Chow test [48] of NDVI and time was used to confirm whether NDVI phases presented structural differences at the selected breaking points.

Correlations of NDVI with Meteorological Variables
The NDVI data were based on eight-day compound images. The image was selected based on criteria such as clouds and solar zenith. The temperature was the average for every eight days. Precipitation was accumulated every eight days, being coherent with the time scale of NDVI series. For exploring the existence of temporal patterns, we focused on two types of correlations: Pearson's correlation coefficient of the NDVI values with Temperature (Temp) and Precipitation (Pp) at different phases (i) is: where V i (t) is the average of the 18 years of the variable V (NDVI, Temp or Pp) at time t belonging to phase i; Pearson's correlation coefficient of the NDVI series belonging to phase i, through the 18 years (j), with each of the climatic variables (Temp and Pp) at different time lags (s) is: where NDVI i (j, t) are the NDVI values at year j and time t that belong to phase i. The Temp(j, t − s) are the temperature values at year j and delayed s times the lag time, which is eight days. Analogously, Pp(j, t − s) can be defined.

Aridity Index and NDVI
The aridity index was calculated following [35], but instead of accumulating annually, we used the phases described by the NDVI patterns: where P i is the summation of the accumulated precipitation of each phase for each year and was analogously done for the accumulated potential evapotranspiration ( ETo i ). Then, the average NDVI value for each phase was calculated. The cumulative aridity index and the cumulative average NDVI for each phase were plotted, and a linear regression was calculated to compare the four areas. A high slope would indicate an efficient use of its water resources.

Interannual Variation
When yearly average NDVI, temperature, and precipitation were plotted at the two resolutions, different behaviors were found in these 18 years (Figures 3 and 4). In all the areas, NDVI at both resolutions were very similar. In A1 and A3, NDVI presents almost a constant value (A1 has 20 and NDVI nearly 30) as does temperature (15 • C in every area, except A3 where it reaches 18 • C). Precipitation is constant in A1 while in A3 it shows an increasing trend (A3 presents the lowest precipitation, mostly below 300 mm and the remaining areas range between 300 and 400 mm). On the other hand, A2 and A4 NDVI present a slight increase (with values around 20 for A2 and 40 for A4). In both areas, temperature and precipitation show different trends.
Remote Sens. 2021, 13, x FOR PEER REVIEW 6 present a slight increase (with values around 20 for A2 and 40 for A4). In both ar temperature and precipitation show different trends.  To study whether these trends, visible in Figures 3 and 4, are statistically signific a Mann-Kendall test [49,50] was applied in each area (Table 1). A1 shows a decrea trend with NDVI MR and LR, but it is only significant with MR. On the other hand, A A4 show an increasing trend, significant in both resolutions. The temperature incre in A1 to A3 and decreases in A4 and precipitation decreases in A1 and A4 and incre in A2 and A3. However, temperature and precipitation do not present significant tre except for precipitation in A4, which decreases significantly (Table 1).

Areas
Trend Remote Sens. 2021, 13, x FOR PEER REVIEW 6 present a slight increase (with values around 20 for A2 and 40 for A4). In both ar temperature and precipitation show different trends.  To study whether these trends, visible in Figures 3 and 4, are statistically signific a Mann-Kendall test [49,50] was applied in each area (Table 1). A1 shows a decrea trend with NDVI MR and LR, but it is only significant with MR. On the other hand, A A4 show an increasing trend, significant in both resolutions. The temperature incre in A1 to A3 and decreases in A4 and precipitation decreases in A1 and A4 and incre in A2 and A3. However, temperature and precipitation do not present significant tre except for precipitation in A4, which decreases significantly (Table 1).

Areas
Trend To study whether these trends, visible in Figures 3 and 4, are statistically significant, a Mann-Kendall test [49,50] was applied in each area (Table 1). A1 shows a decreasing trend with NDVI MR and LR, but it is only significant with MR. On the other hand, A2 to A4 show an increasing trend, significant in both resolutions. The temperature increases in A1 to A3 and decreases in A4 and precipitation decreases in A1 and A4 and increases in A2 and A3. However, temperature and precipitation do not present significant trends, except for precipitation in A4, which decreases significantly (Table 1).

Low and Medium Resolution NDVI Series
The average NDVI time series with LR and MR are shown in the left column of Figures 5 and 6 for the four areas studied. When the averages of all pixels for both resolutions are plotted, they tend to coincide. Nevertheless, we can differentiate them at the peaks and valleys of most years. In A1 and A2 ( Figure 5), the two series are very similar. However, we can see some peaks where LR is higher for A1 and A2. This effect is more prominent in the A3 case. On the other hand, in A4, the differences are minor between the two resolutions. These peaks where LR is higher than MR are due to a few pixels and dates where the NDVI are much higher in both resolutions. However, because the MR has more pixels averaged, the smoothing effect is more noticeable. Since not all pixels were used for A3, this rise of LR is more conspicuous.

Low and Medium Resolution NDVI Series
The average NDVI time series with LR and MR are shown in the left column of Figures 5 and 6 for the four areas studied. When the averages of all pixels for both resolutions are plotted, they tend to coincide. Nevertheless, we can differentiate them at the peaks and valleys of most years. In A1 and A2 ( Figure 5), the two series are very similar. However, we can see some peaks where LR is higher for A1 and A2. This effect is more prominent in the A3 case. On the other hand, in A4, the differences are minor between the two resolutions. These peaks where LR is higher than MR are due to a few pixels and dates where the NDVI are much higher in both resolutions. However, because the MR has more pixels averaged, the smoothing effect is more noticeable. Since not all pixels were used for A3, this rise of LR is more conspicuous.  In these years, there are hotter temperatures for negative values with reduced precipitation. When ZNDVI goes above 2, high peaks on precipitation and temperature are lower than in other years. ZNDVI in A3 only goes above two once, at the end of 2019, coinciding with a strong peak in precipitation (Figures 4 and  6). Area 4 does not show values beyond these levels, although it shows how the trend rises or drops with these events. Regarding the resolution in A1 to A3, we can see peaks and valleys where LR has a higher frequency of extreme values. These differences are smoother for A4.

NDVI Patterns and Meteorological Variables
The LR and MR NDVI annual evolution and temperature and precipitation are shown in the following boxplots for each of the studied areas (Figures 7-10). The LR and MR NDVI values in A1 and A2 are very similar, then A3 presents an increase in these values, and A4 has much higher NDVI values than the rest. This pattern may be, at leas in part, explained by increasing tree coverage in the areas. In each graph, the year was divided into different phases as the NDVI average trend changed. In these years, there are hotter temperatures for negative values with reduced precipitation. When Z NDVI goes above 2, high peaks on precipitation and temperature are lower than in other years. Z NDVI in A3 only goes above two once, at the end of 2019, coinciding with a strong peak in precipitation (Figures 4 and 6). Area 4 does not show values beyond these levels, although it shows how the trend rises or drops with these events. Regarding the resolution in A1 to A3, we can see peaks and valleys where LR has a higher frequency of extreme values. These differences are smoother for A4.

NDVI Patterns and Meteorological Variables
The LR and MR NDVI annual evolution and temperature and precipitation are shown in the following boxplots for each of the studied areas (Figures 7-10). The LR and MR NDVI values in A1 and A2 are very similar, then A3 presents an increase in these values, and A4 has much higher NDVI values than the rest. This pattern may be, at least in part, explained by increasing tree coverage in the areas. In each graph, the year was divided into different phases as the NDVI average trend changed.
A1 and A2 present a larger dispersion during the spring (March, April, and May). This dispersion is present but less prominent in A3, probably due to its less abundant precipitation. A4 exhibits a more consistent dispersion of values throughout the year with a reduction during the summer (June, July, and August), given by the almost complete tree coverage of this area. NDVI values decrease during the summer in all areas with a more visible trend in A4. When comparing NDVI and meteorological variables, the former shows greater heterogeneity among areas. NDVI and temperature appear to show a delay between their peak values. The average temperature peaks occur on 28 July in all the zones. The NDVI peaks on the 20 July in all the areas except A4, which showed it on 12 July (Figures 7 and 9 for A1 and A2, and Figures 8 and 10 for A3 and A4). Area A1 to A3 are more heavily influenced by agricultural practices, which may cause the delay between NDVI and temperature, which only takes eight days. Area A4, as open forest, does not have any irrigation regimes that can relate to an earlier peak of NDVI when temperatures rise. No different trends were found when boxplots were plotted for weeks, fortnight,        A1 and A2 present a larger dispersion during the spring (March, April, and M This dispersion is present but less prominent in A3, probably due to its less abund precipitation. A4 exhibits a more consistent dispersion of values throughout the with a reduction during the summer (June, July, and August), given by the almost c plete tree coverage of this area. NDVI values decrease during the summer in all a with a more visible trend in A4. When comparing NDVI and meteorological varia the former shows greater heterogeneity among areas. NDVI and temperature appea show a delay between their peak values. The average temperature peaks occur on Jul in all the zones. The NDVI peaks on the July 20 in all the areas except A4, which sho it on July 12 (Figures 7 and 9 for A1 and A2, and Figures 8 and 10 for A3 and A4). Are to A3 are more heavily influenced by agricultural practices, which may cause the d between NDVI and temperature, which only takes eight days. Area A4, as open fo does not have any irrigation regimes that can relate to an earlier peak of NDVI w temperatures rise. No different trends were found when boxplots were plotted for we fortnight, months, or season (data not reported). Meteorological trends remain sim Because the NDVI trends did not match the change of climatic seasons, we separated the NDVI and meteorological variables according to the NDVI trend changes, as shown in Figures 7-10. Areas 1-3 were split into five phases, whereas A4 was split into four phases ( Table 2).

Intra-Annual Regression by Phases
The regression analysis results between NDVI with temperature and precipitation by phases using the average of 18 years are shown in Figures Table 2). The regression analysis shows that for A1 to A3 (Figures A1-A3, Appendix B) temperature has a more substantial effect in NDVI in phase 2, 3, and 5, compared to phases 1 and 4, although phase 5 shows a slightly weaker relationship for A2, and especially A3. In A4, phases 3 and 4 show a strong relationship in the regression analysis and mild for phase 2 ( Figure A4 in Appendix B).
The regression analysis for precipitation shows that for A1 to A3 ( Figures A5-A7 in Appendix B), the herbaceous areas, NDVI in phase 3 is influenced by precipitation, while the others show a fragile relationship. No strong relationships are found in any phases in A4 ( Figure A8). Due to the use of average values in the regression, differences between the two resolutions are small, although LR tends to have lower R 2 values than MR.

NDVI and Meteorological Series Correlation
The Chow test confirmed significant structural differences between phases for all areas (Table 3). Medium-and low-resolution NDVI produced similar results in the previous analyses. However, there were differences, particularly regarding the lagged responses. This paragraph describes the results of MR NDVI correlation analysis and discusses the differences with LR NDVI. The correlation values mentioned are the highest lagged correlation (Tables 4 and 5). NDVI has similar lags in each phase throughout A1 to A3. In phases 1 and 4, the NDVI values experience little change (with a difference within these phases between 4 and 10 in NDVI).   For this reason, temperature and precipitation have a minimal correlation with NDVI during these phases. Phase 2 shows a positive correlation of NDVI and temperature for these three areas: 0.32 for A1, 0.28 for A2, and 0.25 for A3. Phase 3 presents a moderate negative correlation for A1 and A2 (−0.73 and −0.61, respectively) and a low correlation (−0.12) in A3. Phase 5 has a low negative correlation for A1 (−0.32), A2 (−0.35), and A3 (−0.04). Area 4 always presents a negative correlation between temperature and NDVI, all lower than −0.4 except in phase 3, which has a stronger negative correlation of −0.73. Precipitation presents small positive values in all phases and areas, all below 0.35.
Lags on temperature for A1 to A3 in phase 1 and 4 range from 0 to 8 days, except for phase 4 in A2. All of these have very low values due to the constant behavior of NDVI values. In phase 2, the highest correlation between NDVI and temperature is found with a 16-day lag in A1, and 32-day lag in A2 and A3. Phase 3 has lags of 24 and 32 days for these three areas. Phase 5 differs more between areas, with a lag of 32 days for A1, eight days for A2, and 24 days in A3. Area 4 only shows a 16-day lag in phase 3, while the Remote Sens. 2021, 13, 840 13 of 24 rest of its phases present no lag. The differences between the areas may be related to two factors, differences between their vegetation and their hydrological regimes. Precipitation correlation lags show small similarities between all four areas. The scattered precipitation of these habitats and water scarcity makes it difficult to find a pattern across the areas. These disparities appear to be related to differences in heavy rains at specific times in each phase and area.
Comparing these results with the LR ones (Tables 6 and 7), correlation exhibited one of our analyses' significant differences between the two resolutions. The lag presented by the highest correlation with temperature was similar for MR and LR, except for some cases where there was an approximate difference of 8-days lag. However, phases 4 and 5 present more significant differences with a lag of 16 and 32 days, in some cases. The precipitation variable presents an almost constant smaller lag in the LR, ranging from 8 to 48 days. Partial correlation with NDVI, temperature, and precipitation was made, but no significant differences were found (Tables A2-A5, in Appendix C). They indicated that water availability is highly affected by temperature in this region.   Figure 11 shows the four areas cluster into two groups. On one side, we find A1 and A2 intermingling. These two areas present lower cumulative NDVI than A3 and A4, and higher cumulative values of AI than A4 and A3. The AI for A4 was calculated based on its four phases, instead of five as was done for the other areas. However, this change does not affect the pattern and slope of this graphic. In the upper part of Figure 11, we find A3 and A4 showing efficient use of water resources. Their higher slope reflects the increase of NDVI as AI rises. With more typical Mediterranean vegetation, A3 and A4 have more efficient water use, particularly A3 with xerophilous shrubs as opposed to less xeric shrub vegetation, such as rosemary found in A4.  Figure 11 shows the four areas cluster into two groups. On one side, we find A1 and A2 intermingling. These two areas present lower cumulative NDVI than A3 and A4, and higher cumulative values of AI than A4 and A3. The AI for A4 was calculated based on its four phases, instead of five as was done for the other areas. However, this change does not affect the pattern and slope of this graphic. In the upper part of Figure 11, we find A3 and A4 showing efficient use of water resources. Their higher slope reflects the increase of NDVI as AI rises. With more typical Mediterranean vegetation, A3 and A4 have more efficient water use, particularly A3 with xerophilous shrubs as opposed to less xeric shrub vegetation, such as rosemary found in A4.

NDVI Series Comparison
The average NDVI values for LR and MR of MODIS display relatively strong simi larities. However, differences between them on their series' peaks and valleys can be detected (Figures 5 and 6). These differential patterns among areas suggest that it is caused by the difference of the averaged pixels within each area, showing that a finer scale is more representative of the area. On the other hand, some peaks and valleys are more pronounced, mainly when extreme meteorological events occur. These differences show a spatial difference among the pixels, which should be further researched.
These differences are confirmed when we compared the boxplots and the results o Pearson correlation and regression analysis. We found the significant differences in the boxplots among resolutions in A4, which is the most heterogeneous landscape, with a mix of tree-, scrub-, and grassy-dominated pixels. The medium resolution provided more significant results when studying their tendencies; the Mann-Kendall test for A1 with MR showed a tendency that was not visible with LR. Regression analysis provided more

NDVI Series Comparison
The average NDVI values for LR and MR of MODIS display relatively strong similarities. However, differences between them on their series' peaks and valleys can be detected (Figures 5 and 6). These differential patterns among areas suggest that it is caused by the difference of the averaged pixels within each area, showing that a finer scale is more representative of the area. On the other hand, some peaks and valleys are more pronounced, mainly when extreme meteorological events occur. These differences show a spatial difference among the pixels, which should be further researched.
These differences are confirmed when we compared the boxplots and the results of Pearson correlation and regression analysis. We found the significant differences in the boxplots among resolutions in A4, which is the most heterogeneous landscape, with a mix of tree-, scrub-, and grassy-dominated pixels. The medium resolution provided more significant results when studying their tendencies; the Mann-Kendall test for A1 with MR showed a tendency that was not visible with LR. Regression analysis provided more robust results for the relationship of NDVI and climatic variables with MR in A2 and A3, while A1 and A4 had smoother results than A2 and A3. Additionally, MR NDVI had more delayed lagged correlations. The use of medium or low resolution may depend on the spatial heterogeneity of the areas of study. This is probably due to a more detailed depiction of NDVI values in MR, as a smaller pixel size allows us to separate areas that could be spatially diverse. MR allows emerging a more lagged correlation without being tampered by the surroundings if an area dries at different times.
Using both resolutions shows a clearer image of the selected areas, especially highlighting those that are more heterogeneous. These results are in agreement with [31]. The use of Z NDVI and its comparison to the NDVI series also highlights the more significant changes during the series and the differences between the two resolutions. In particular, we find that extreme events are more pronounced when using LR as compared to MR.

NDVI Data and Meteorological Variables
NDVI values are highly related to water availability, as stated by [51]. NDVI responses are strongly linked to temperature in Mediterranean habitats, although this relationship weakens when precipitations are high [52]. In Murcia, the precipitations are scarce, with dry winters and almost no precipitation during the summer months. Furthermore, the temperature rises to its peak in July and August, leading to decreasing NDVI values. These values rise when the precipitation begins, and the temperature starts to descend. Our data show similar trends of higher NDVI values when precipitation increases, as highlighted by [53,54]. Area 4, located in the northwest county of Murcia, has relatively more water abundance and larger NDVI values.
Our correlation analysis shows that NDVI data are strongly and negatively correlated to average temperature when precipitation is strongly limited. The intraseasonal variation of the relationship of the NDVI and the meteorological variables has been approached by using different phases based on the NDVI behavior. The patterns of NDVI correlations with temperature change depending on the land use and the phase selected. Phases 1 and 4 show weak correlation values because the NDVI presents a steady pattern. On the other hand, phases 2, 3, and 5 have either increasing or decreasing NDVI trends. These phases present negative values except for phase 2 for A1 to A3 when the spring precipitation occurs. The correlations for NDVI and precipitation are very low, all below 0.35. However, the regression analysis highlights a stronger relationship for A1 to A3 for phase 3, when the precipitation starts high and steadily drops as temperature increases. A4, covered with trees, does not show a strong relationship between NDVI and precipitation where its precipitation in phase 2 is more limited than in the other areas. These changes among phases suggest that one of the critical factors affecting NDVI is water availability, matching other areas with semiarid and arid climates [55,56].
After characterizing the climatic variables and NDVI dynamics, we approached how the NDVI was affected by the AI, as a combination of precipitation and potential evapotranspiration, reflecting the temperature as well. This relation was established by accumulating NDVI and AI by the phases defined in each area. This allowed us to see how the different vegetation types representing each area respond differently to water availability. A3 and A4 show a more efficient response of NDVI to the rise of AI. Among these two areas, A3 exhibits a better response, in agreement with its more xeric vegetation, particularly in the selected top half pixels. No trees are found and xerophilous scrubs, such as esparto, inhabit the area. On the other hand, A1 and A2, representing different types of crops, show a smaller NDVI in response to a similar AI than A4, and higher than A3.

Conclusions
Comparison of the two scales of MODIS MOD09Q1.006 (MR) and MOD09A1.006 (LR) show that NDVI is scale dependent. These resolutions show differences, particularly when studying correlation and regression analysis. They are suggesting that medium resolution is more suited for spatial and lagged temporal patterns. However, when averaged, the trends are similar between these two resolutions. Lower resolution scales can be used when the studied areas are not spatially heterogeneous for temporal trend analysis, but larger resolutions scales are recommended on spatially diverse areas.
Among the climatic variables used, temperature shows the strongest relationship with average NDVI. Our results reveal that complex interactions of precipitation and temperature may explain real-time NDVI evolution. However, their behaviors vary across the phases. The use of phases based on NDVI patterns, instead of seasons, allows us to describe a more realistic depiction of the arid environments, based on their vegetation dynamics. The study shows a strong positive correlation of NDVI with temperature when high precipitation occurs. Precipitation, however, shows a weak correlation with NDVI. The behavior of both climatic variables points to water availability as one of the major drivers of NDVI in Murcia, as it is suggested by the positive correlations of temperature and NDVI during phases where heavier precipitation occurs.
Aridity index and NDVI allowed us to cluster our four areas into two large groups, A1 and A2; mainly grazed wastelands showed a low increase of NDVI with a high aridity index. In contrast, A4 and A3 with a similar or lower aridity index presented a higher NDVI accumulation, showing more efficient use of available water; A3, especially, presented a higher slope than A4. However, more rangelands and other ecosystems should be analyzed to determine whether or not these differences can discriminate and characterize among other types of rangelands.
Intraseasonal relationship of NDVI with climatic variables was studied by splitting the analyses according to NDVI patterns. This allowed for viewing NDVI dynamics that were obscured when the seasonal division was followed. Intraseasonal and interseasonal characteristics should be taken into consideration in the definition of agrometeorological indexes in rangelands. This study provides a discriminating technique for rangeland management and policymakers. It is expected that future research will expand the knowledge of NDVI drivers at different scales, to develop tools and indexes that can help further comprehend vegetation communities of agricultural lands.  Acknowledgments: The data provided by Appears Team and AEMET are greatly appreciated.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.               Lag Time(Days) Figure A8. Regression analysis of average NDVI and precipitation series separated by phases for LR (a) and MR (b) for A4.