Climate and aridity measures relationships with spectral vegetation indices across desert fringe shrublands in the South-Eastern Mediterranean Basin

Mediterranean regions are hot spots of climate change, where the expected decrease in water resources threatens the sustainability of shrublands at their arid margins. Studying spectral vegetation indices’ relationships with rainfall and potential evapotranspiration (PET) changes across Mediterranean to arid transition zones is instrumental for developing methods for mapping and monitoring the effects of climate change on desert fringe shrublands. Here we examined relationships between 17 spectral vegetation indices (VIs) and four climate and aridity measures, i.e., rainfall, PET, aridity index (AI), and water deficit (WD), calculated at accumulation lags between 1 and 6 months. For this purpose, VIs for 38 sites (100 × 100 m each) representing less disturbed areas were extracted from Sentinel 2A images for 3 years with high (2016), low (2017), and average (2018) annual rainfall. Most of the VIs had shown the highest correlation with the four climate and aridity measures at 2-month accumulation interval. While NDVI relationships with climate measures gained the widest use, our data suggest that indices combining NIR and SWIR bands better correlate with climate parameters. AI is one of the leading annual measures of dryness worldwide; when calculating it monthly, WD was found to represent better the balance between precipitation and PET across the climate transition zone and to be better correlated with VIs. Relationships between NIR and SWIR VIs and water deficit may thus facilitate improvements in monitoring and mapping desert fringe shrublands’ responses to climate change if supported by similar results from other areas.


Introduction
Under current global climate change trends, extensive semi-arid and desert fringe ecosystems experience declining water resources and increasing atmospheric temperatures (e.g., Blondel et al., 2010;Maestre et al., 2012;Polade et al., 2017;Pugnaire, et al., 2019;Tuel & Eltahir, 2020). Since it is introduced by Martonne (1926), the aridity index (AI) with its Thornthwaite (Budyko, 1974;Huschke, 1959) and UNEP (Barrow, 1992;Greve et al., 2019) versions has been widely used to monitor the spread of aridity across these ecosystems (e.g., Pellicone et al., 2019;Speich, 2019;Traore et al., 2020;Derdous et al., 2021;Zomer et al., 2022). Other aridity and dryness indicators were calculated based on the difference between precipitation and potential evapotranspiration with Bussay et al. (2012) and Vicente-Serrano et al. (2012)  as water deficit (WD), while Stephenson (1990), Saha et al. (2020), Abatzoglou et al. (2018), and Brooks (2004) refer to it as climate water balance (CWB). SPEI (Standardized Precipitation Evapotranspiration Index) representing the distribution of deviations from the average water deficit was implemented over broad regions and globally (Vicente-Serrano et al., 2013). Linking spatial and temporal variations of rainfall, PET, AI and WD to surface (soil and vegetation) information derived from satellite images may contribute to a better understanding of the potential environmental impacts from the decline of water resources. Furthermore, statistical formulation of relationships between climatic parameters and imagery spectral properties may serve densifying spatial and temporal information of rainfall, PET, and aridity.
The usefulness of SWIR bands for estimating vegetation water content was suggested first by Tucker (1980) and then demonstrated by Elvidge's (1990) measurements. Gao (1996) claimed that ratios combining NIR and SWIR bands "… failed to detect plant water status change within a biological meaningful range" and explained it by the broad width of Landsat TM SWIR band 5 (1.55-1.75 mm). The better resolution of the SWIR bands of the MODIS sensor yielded successful use of indices combining NIR and SWIR bands, such as normalized difference snow index (Avtar et al., 2020), land surface water index (Chandrasekar et al., 2010), and normalized difference water index (Gu et al., 2007). However, the spatial resolution of MODIS in the SWIR is 500 m limiting the application of these indices to regions with relatively homogenous vegetation cover. High spatial heterogeneity of Mediterranean landscapes due to geological, lithological, topographic, climatic, and anthropogenic diversity (Benito-Calvo et al., 2009;De Falco et al., 2021;Ferrando et al., 2021;Mooney, 1988;Safford & Vallejo, 2019) results in changes of ecological patterns within short distances. Life form and species compositions with different modes of response to rainfall and atmospheric heat would highly vary within pixels of 500 × 500 m (e.g., MODIS). SENTINEL-2 sensor offers improved SWIR spectral resolution compared with the Landsat 5 (0.091 versus 0.2 mm) combined with a good spatial resolution which limits the within-pixel variation of life forms and species compositions on the one hand and increases the differentiation between types of ecological patterns on the other hand.
While most of the research on climatic parameters' relationships with VIs is based on annual average data, a significant number of studies assessed different cumulative/averaging intervals between 1 and 6 months, primarily for rainfall (Al-Bakri & Suleiman, 2004;Chamaille-Jammes et al., 2006;Chen et al., 2014;Daba & Oromia, 2016;Dagnachew et al., 2020;Lin et al., 2016;Nicholson & Farrar, 1994;Nicholson et al., 1990;Verma et al., 2022;Wu et al., 2015). Fewer studies also included cumulative/average values for other climatic parameters: seasonal land surface temperature (Dyosi et al., 2021), atmospheric temperature (Ebrahimi-Khusfi et al., 2020), air temperature, and sunshine duration within the context of climate warming anomalies (Zeng & Yang, 2008), and SPI and SPIE (Mukwada & Manatsa, 2018). The study of relationships between multiple Vis at the VIS-NIR-SWIR range with multiple climatic parameters at different time intervals across desert fringe transitions zones is instrumental for better understanding their plant patterns and potential response to future global warming: it is a "space for time" reasoning commonly used for this purpose (Roitberg & Shoshany, 2017). Hence, the following research questions are addressed: (1) What are the accumulation intervals of four measures (rainfall, PET, WD, and AI) which best correlate the different VIs? (2) Which vegetation index best correlates these four climatic parameters when accumulated within the best interval?
These questions are assessed using VI data for the climatic gradient in central Israel extracted from Sentinel-2A images with a spatial resolution of 10 m/ pixel and 20 m/pixel, for 3 years representing low (2017), average (2018), and high (2016) rainfall.  (Nicholson et al., 1990

Study area
Dry summers and wet winters characterize the five Mediterranean-type climate (MTC) regions, which include ecosystems of the margins of the deserts bordering them (Rundel, 1998;Rundel et al., 2018). The diversity of rainfall, radiation, topography, lithology, and human effects in general, and fires in particular, resulted in high species richness, compared with other regions across the globe. Although the five MTC regions are located on different continents, their floristic and life-form compositions are remarkably similar (e.g., Naveh, 1967;Di Castri, 1981;Shmida, 1981;Carmel & Flather, 2004;Naveh, 2019).
The study area stretches along the rainfall gradient between 500 and 250 mm/year over Rendzina soil developed on chalk of the middle Eocene age with gentle hilly topography (Dan et al., 1976). Sites were located on less disturbed areas based on visual assessment of plant, soil, and rock patterns as observed from 2021 Survey of Israel (SOI) orthophoto (www. govmap. gov. il) in comparison to patterns found in earlier surveys using orthophoto from 2010 (Shoshany & Karnibad, 2015), ultralight flight photographs (Shoshany, 2012) and high-resolution hyperspectral Fenix strip from 2017 (Paz-Kagan et al., 2021). Overall, there were delineated 38 plot-areas sized 100 × 100 m (Fig. 1).

Climate and remote sensing data acquisition and analysis
Our research methodology is comprised of three parts: the first consists of mapping the rainfall and PET across the climatic gradient, facilitating derivation of spatial distributions of WD and AI, and allowing calculation of cumulative values for 1 to 6 months for the four climate/ aridity measures; the second concerns calculating vegetation indices from the Sentinel images, and the third involves an analysis of the relationships between those cumulative values and the spectral vegetation indices.
Rainfall and temperature data were downloaded for this period from the online archives of the Israel Meteorological Service (https:// ims. data. gov. il/ ims/1). Monthly rainfall data were obtained in millimetres/month for 31 meteorological stations, and three hourly temperature readings (in degrees Celsius) were obtained from 27 meteorological stations, and further aggregated into daily and then mean monthly temperatures. Monthly PET values (mm/month) were calculated based on Thornthwaite's method (1948) using the following equation (Xu & Singh, 2001): where c = 16 is a constant, d is the average day length in a given month, N is the number of days in a given 10T I month, T a is the average monthly temperature (C°), I is the annual heat index obtained from Eq. (2), and α is the exponent calculated from Eq. (3).
where i is the number of months (from 1 to 12). Figure 2 presents the monthly rainfall and PET fluctuations along the climatic gradient, expressing the inverse temporal patterns of rainfall and PET, respectively, between the short rainy winters and the long and hot summers.
Continuous coverage of the monthly rainfall and PET from the meteorological stations was obtained using local polynomial (LP) interpolation, provided by ArcMap 10.7.1 software, following its comparison with inverse distance weighting (IDW), original kriging (OK), and universal kriging (UK).
While IDW and UK showed low accuracies, OK with exponential variogram yielded higher accuracies than LP during months with moderate and high rainfall and PET levels. However, OK yielded spatial discontinuities and had lower accuracies than LP in months when rainfall or PET dropped below moderate levels. LP formed consistently continuous rainfall and PET surfaces during all seasons and therefore was preferred for this study. The interpolated grid was resampled with a 10-m/pixel resolution over the entire study area to allow matching with the vegetation indices grid created from the Sentinel-2A imagery.
As discussed earlier, relationships between rainfall and PET attracted considerable attention, especially since the introduction of the AI as the ratio between rainfall and PET and its use in the preparation of the World Atlas of Desertification (UNEP, 1992). The difference between rainfall and PET was used in a lower number of studies to estimate water deficit (e.g., Bussay et al., 2012;Vicente-Serrano et al., 2012) or the CWB (e.g., Abatzoglou et al., 2018;Brooks, 2004;Saha et al., 2020) and was also mapped over wide regions.
(https:// www. usgs. gov/ media/ images/ map-gridd edvalues-1971-2000-avg-preci pitat ion-minus-avg-pet). Cumulative rainfall, PET, WD, and AI were calculated for time spans between 1 and 6 months with a time lag of 1 month prior to the date of the satellite image, using the following equations: where CR ij , CP ij , CWD ij , and CAI ij are cumulative rainfall, PET, and water deficit, and aridity index, respectively, calculated for month i , which corresponds to the month of satellite image capturing, and j is the number of cumulative months (from 1 to 6) preceding the month i . For example, cumulative rainfall of 1 month for a satellite image representing vegetation in December 2015 equals monthly rainfall in November 2015; cumulative rainfall of 2 months for a satellite image representing vegetation in December 2015 equals the sum of monthly rainfall in November and October 2015, and so on.
AI is measured primarily annually, while here it is calculated at different monthly accumulation lags. At the seasonal resolution, WD better represents the seasonal balance between precipitation and PET than AI, as the second measure converges to 0 during summer month.

Calculating vegetation indices
Time series of 25 Sentinel-2A images with spatial resolutions of 10 and 20 m/pixel taken between December 2015 and November 2018 were acquired for this study from https:// scihub. coper nicus. eu/ dhus/#/ home. All images were downloaded at Level-1C, and then atmospherically corrected using the DOS1 algorithm from the semi-automatic classification plugin in QGIS. Only images with cloud cover lower than 10% were acquired ( Table 2). The 20-m and 20-m spatial resolution bands were resampled at 10 m resolution, and 17 vegetation indices were calculated. Table 3 presents the equations for each of these indices and Fig. 3 presents monthly variations in four representative vegetation indices for four sites.

Analysis of vegetation indices relationships with climatic measures
Linear regression equations were calculated between each site (100 × 100 m) average VIs and its corresponding cumulative values of rainfall (CR), PET (CPET), WD (CWD), and AI (CAI) for six accumulation intervals (between 1 and 6 months) for each month (m) between December 2015 and November 2018.
Following preliminary assessment, relationships between CPET and vegetation indices were found to be non-linear with best linear fit to the log-normal of CPET: where s is the site index (1 through 38). Figure 4 presents the data processing flow where in the inner core loop there is calculated one of the regression Eqs.(8)-(11) and its coefficient of determination (R 2 ), for each site. The scatter diagrams of average R 2 for all sites and their standard deviation as determined for each of the VI's are presented on the left side of Fig. 5 for each of the climate/aridity measures. The statistical relationship between a vegetation index and a climatic measure is stronger when R 2 is high and when the standard deviation of its variation is low. Assessment of model performance by combining correlation and standard deviation was proposed by Taylor (2001); however, here we look at the statistics of a distribution of R 2 rather than the correlation between observations. Thus, a new combined parameter, square root of average R 2 and 1-Std (SRS), was formalized as follows: SRS values range between 0 and 1.414 (

Results
Of the 17 indices, 10 for the rainfall, 14 for the PET, and 14 indices for WD were having the highest SRS for a 2-month accumulation interval (Table 4, Fig. 6). Highest SRSs were found on the other hand for a 3-month accumulation interval in 10 indices for AI. However, the difference between the 2-month and 3-month SRS averages was only 0. 01 (1.10 for a 3-month versus 1.09 for a 2-month).
Ranking the vegetation indices according to their SRS with cumulative rainfall, PET, AI, and WD for 2 months (Fig. 7) indicates that the same seven vegetation indices (GARI, MSI1610, MSI2190, NDII1610, NDII2190, NDVI, and RENDVI) had the highest SRS values.
Four NIR-SWIR indices, MSI 1610 , MSI 2190 , NDII 1610 , and NDII 2190 together with NDVI were consistently ranked highest in terms of SRS for all four climatic parameters. In between the four climate/aridity measures, WD has the highest SRS values and AI has the lowest, with all indices below SRS of 1.2 (Table 4 and Fig. 7).
The slope and intercept coefficients of the linear regressions calculated between the VIs and climatic parameters (Eqs. (8)-(11)) for each site represent the spectral reflectance response to changes in the sites' climatic and habitat conditions and their species compositions along the gradient. VIs with the widest variation in slope and intercept between sites thus best represent, hypothetically, ecological differences across the climatic gradient. Figure 8 (a) shows that the scatter of slope and intercept values for the best (highest SRS) NIR-SWIR VI (MSI1610) is wider than that of the corresponding scatter for the VIS-NIR VI with the highest SRS (NDVI). Calculating the standard deviation for the coefficients of the seven leading VIs (Fig. 8b) indicates that MSI1610 and MSI2190 exhibit significantly higher variation than all other indices and, thus, may be expected to better serve in the mapping of vegetation conditions along the semi-arid to an arid climatic gradient.  Huete et al. (2002) Green atmospherically resistant index Gitelson et al. (1996) Green norma. difference vegetation index GNDVI = NIR−GREEN NIR+GREEN Gitelson and Merzlyak (1998) Global vegetation moisture index GVMI = (NIR+0.1)−(SWIR 1610nm +0.2) (NIR+0.1)+(SWIR 1610nm +0.2) Ceccato et al. (2002) Moisture stress index at 1610 nm MSI 1610nm = SWIR 1610nm NIR Hunt and Rock (1989) Moisture stress index at 2190 nm MSI 2190nm = SWIR 2190nm NIR Hunt and Rock (1989) Normalized diff. infrared index at 1610 nm

NIR+[GREEN− (BLUE−RED)]
Hardisky et al. (1983) Normalized difference infrared index at 2190 nm Hardisky et al. (1983) Normalized difference vegetation index NDVI = NIR−RED NIR+RED Rouse et al. (1973) Normalized multi band drought index Wang and Qu (2007) Soil adjusted vegetation index OSAVI = NIR−RED NIR+RED+0.16 Rondeaux et al (1996) Red edge inflection point   African savanna, by Chen et al. (2014) for mainland Australia, by Wu et al. (2015) for middle-and low-latitude shrubs, and by Daba and Oromia (2016) for southern Ethiopia. Two-month accumulation intervals were suggested to be best by Nicholson et al. (1990) for the Sahel and East Africa, by Nicholson and Farrar (1994) for Botswana, and by Lin et al. (2016) for the Liao River Basin in Jilin Province, China. Dagnachew et al. (2020) found that the highest correlation between NDVI and rainfall is obtained for rainfall lagged by 2 or 3 months in the Gojeb River Catchment, Omo-Gibe Basin, Ethiopia. Durante et al. (2009) suggested that land use determines the best accumulation times, reaching up to 4 months for Mediterranean Iberian ecosystems. Vicente-Serrano et al.
(2013) found that for dry forests, shrublands, and steppes, optimal accumulation times are between 2 and 4 months, while for semi-arid and semi-humid biomes there is a need for longer time scales (8-10 months). Studies that assessed temperature and evapotranspiration in addition to rainfall with accumulation lags equal to and longer than 3 months were reported by Chuai et al. (2013) for inner Mongolia, China, and by Gouveia et al. (2017) for the Mediterranean basin.
To improve our understanding of VI responses to accumulated climatic conditions across transition zones between semi-arid and arid regions, we expanded the research on the affinity of climatic parameters with spectral vegetation indices, by including PET, AI, and WD While rainfall is limited to the short rainy season and converges in all sites to almost zero during the summer, PET temporal variation is highest during the long dry season and is non-zero during the winter (Fig. 2). Restricting the analysis of VI's relationships with climatic measures to rainfall thus limits the representation of the climatic conditions determining growth and survival of shrubs in semi-arid to arid transition zones. Utilizing cumulative values compensates to a certain extent the representation of the effect of summer conditions on desert fringe plants. The higher correlation of PET with the VI's compared with rainfall ( Fig. 6) can be thus explained by the additive information content provided by the PET for the spring, summer, and autumn seasons. Water deficit further improves the representation of seasonal spectral plant conditions by combining monthly rainfall and PET information. Such improvement is expressed by higher correlations of the VI's with WD. The aridity index, although it combines rainfall and PET, converges to almost zero during the summer season and thus does not represent the decline in plant's water resources during the hot season and its margins and consequently yield the lowest correlation levels with the different VIs.
The study area between the Judean mountains and the Negev Desert, as other semi-arid and desert fringe zones across Mediterranean regions, underwent numerous human disturbance and recovery cycles during the last centuries (Rundel, 1998;Rundel et al., 2018), resulting in high ecological and geomorphic divergence. Our research sites were selected to represent areas less affected by fires, grazing, woodcutting, and historical agricultural land uses. Yet there remains between-site variation in plant conditions due to local disturbance history and local differences in lithology, soil, topography, and height above the regional water table. All of these sources of surface conditions partly explain the between-site and monthly variation in spectral vegetation indices which deviate from the generalized regression line representing the change in climatic conditions. The result according to which a 2-month accumulation lag yields high R 2 and SRS levels for most of the indices and most of the sites suggests that most of the between-seasons and between-years variations in rainfall, PET, AI, and WD may be well generalized by the locally fitted linear regression models within this accumulation interval. Changes in coefficients of these linear regression models can be expected to express the determinants of the spatial habitat variations across our desert fringe shrublands and steppe.
In general, the regression R 2 and SRS scores for the seven top-ranking indices were relatively high and similar for rainfall, PET, AI, and WD, with the highest correlation for the WD and lowest for the AI. The four VIs representing combinations of NIR and SWIR are included in this group of high-ranking indices, and two or three of them are also among the best four indices for all three climatic parameters. There are multiple environmental, habitat, and local factors governing water utilization by plants in water-limited ecosystems (e.g., Dralle et al., 2020), consequently determining green growth (cover and biomass). The consistent high ranking of NDVI is related to its sensitivity to the relative vegetation cover and the photosynthetic activity. Indices combining NIR and SWIR bands were found to have a higher sensitivity to water content both in plants and soil than NDVI (e.g., Ceccato et al., 2001Ceccato et al., , 2002Tian & Philpot, 2015), a fact which might explain their higher correlation with the four climate/aridity measures. The highest MSI1610 and MSI2190 indices' regression slope and intercept variability across the desert fringe presents to our opinion their sensitivity to net water availability from rainfall infiltration and water accumulation in the root zone following the loss of water due to runoff and evaporation from the atmosphere and the soil. Such sensitivity to variability in water resources may be instrumental for monitoring spatial and temporal changes in desert fringe environments due to climate change and or human disturbance.

Conclusions
The climatic gradient between the Judean Lowland and the Negev Desert represents a transition zone between dense Mediterranean shrublands, open shrublands, desert fringe Batha, and arid ecosystems, which characterizes wide desert margins in North Africa, the southern and eastern Mediterranean Basin, Western Australia and Chile, and California. Understanding the expected response of these water-limited ecosystems to climate change presents a significant and important challenge (e.g., Gordo & Sanz, 2010;Gouveia et al., 2017;Malhi et al., 2020). Studying relationships between vegetation patterns and rainfall and PET is fundamental for this purpose, and remote sensing using spectral vegetation indices is instrumental for investigating these relationships over large regions. Examination of NDVI and rainfall relationships appears to be the main approach, implemented over desert and semi-arid regions (Nicholson et al., 1990;Nicholson & Farrar, 1994;Al-Bakri & Suleman, 2004;Chamaille-Jammes et al., 2006;Chamaille-Jammes & Fritz, 2009;Svoray & Karnieli, 2011;Dagnachew et al., 2020). In this study, we aimed to broaden the examination of these relationships beyond NDVI by using many other VIs and specifically those based on combinations of NIR and SWIR bands. The inclusion of potential evapotranspiration for this purpose was reported in a most limited number of previous studies (Dorman et al., 2013;Gouveia et al., 2017;Islam & Mamun, 2015;Lamchin et al., 2018;Munson et al., 2016;Thoma et al., 2016;Vicente-Serrano et al., 2013). Few studies presented results regarding the time-lag response of VIs to precipitation and PET and the balance between them (AI and WD).
First, this study reveals that for most of the indices and for most of the sites, the seasonal, and annual cycles of rainfall, PET, AI, and WD relationships with corresponding spectral VIs are well described by the locally fitted linear regression models for the 2-month accumulation interval. The second most meaningful result concerns the better statistical relationships of indices combining NIR and SWIR with climate/aridity measures compared with indices combining VIS and NIR. The third important finding concerns the differences between the climate/aridity measures (AI and WD) affinity with desert fringe vegetation conditions as expressed by their correlations with the different VIs. The highest correlation of WD with MSI1610 and MSI2190 compared to other combinations of climate/aridity measures with other VIs and the highest variation in the linear regression coefficients of WD with these two indices suggest that their relationships are highly sensitive to spatial desert fringe shrublands spatio-temporal variation and their potential response to climate change.