The present study area is located in the segment of the arid western plain zone (ACZ-II) of western Rajasthan, India. The geographical extension of the study region is between 74°3′6.316″E and 74°39′9.247″E and 27°47′53.022″N and 28°40′21.049″N encompassing an area of 4,500 km² (Figure 1). There are four main agro-climatic zones in western Rajasthan: (1) irrigated north-western plain (ACZ-I), (2) arid western plain (ACZ-II), (3) transitional plain of inland drainage (ACZ-III), and (4) transitional plain of Luni Basin (ACZ-IV). A negligible fraction of the arid western zone is covered by vegetation (Kar, 2014). The extreme range of temperature (–1°C–49°C) and low rainfall (350 mm) are some of the significant climatic characteristics of this region (Census of India, 2011).

In the present study, 21 years of mean monthly continuous data (2000–2021) for NDVI, precipitation (mm/day), soil moisture (kg m^–2), evapotranspiration (kg m^–2), and 2-m air temperature (^oC) were downloaded from NASA-supported Giovanni datasets. A NASA Goddard Earth Science Data and Information Services Center (GES DISC) Distributed Active Archive Center (DISC) web application (https://giovanni.gsfc.nasa.gov/giovanni/) was used for accessing, visualizing, and analyzing considerable volume geospatial data. The details of the data are shown in Table 1.

Karl Pearson’s coefficient of correlation “r” is a widely used technique for determining correlation among variables (Zeren Cetin et al., 2023). In the present study, multiple correlations between the variables were established and are portrayed in Figure 2 using the following expression (1) in RStudio using the metan package (Olivoto and Lúcio, 2020): r x y = ∑ i = 1 n x i − X ′ y i − Y ′ ∑ i = 1 n x i − X ′ 2 ∑ i = 1 n y i − Y ′ 2 , (1)

where n is the total number of samples, X′ and Y’ are the respective means of x and y, and r_xy is the correlation coefficient between variables x and y.

To determine the trend in the time series of the selected variables, Mann–Kendall (MK) statistics and Sen’s slope were performed in RStudio (Team R. C., 2020) using the Kendall (McLeod, 2016) and trend (Pohlert, 2023) packages. The MK test is a widely used tool for the depiction of monotonic trends in time-series data (2 and 3): S = ∑ i − 1 n − 1 ∑ j = i + 1 n s g n y j − y i , (2)

where s g n y j − y i = s g n R j − R i = f y = − 1 , i f   y j < y i , 0 , i f   y j = y i , + 1 , i f   y j > y i , (3) where S is the MK statistics or Kendall’s tau (MK-tau) and R_j and R_i are the ranks of time series observation of y_j and y_i, respectively. The test assumes that data are independent and identically distributed. The following expressions (4 and 5) can be used for the denervation of mean E(s) and variance var(s): E s = 0 , (4) Var s = n n − 1 2 n + 5 − ∑ p = 1 q t p t p − 1 2 t p + 5 18 , (5)

where t_p is the number of ties for the pth extent and q is the number of tied groups. The second part represents the tuning for tied observations. It is assumed that with a significant number of observations, S tends to follow a normal distribution. The Z significance test (Z_MK) with a p-value <0.05 was computed by Eq. 6: Z M K = S − 1 V a r   S , i f   S > 0 , 0 , i f   S = 0 , S + 1 V a r   S , if  S < 0 . (6)

Here, the null hypothesis (H_o), that is, no trend, is tested against the alternative hypothesis (H₁), having a trend in the time-series data. In addition, Sen’s slope was estimated using the following expression: Q i = X j − X k j − k f o r   i = 1 … … … N , (7)

where Q_i is the estimated slope for each pair of observations; j and k denote time steps where j>k.

The time-series data for all the variables were represented in Origin software. The collected data were categorized in two ways. The entire dataset was initially organized according to the seasons (winter, pre-monsoon, monsoon, and post-monsoon) and represented through a box plot in RStudio using the tidyverse package (Wickham et al., 2019). In addition, the entire dataset was organized in mean monthly results and represented using a boxplot. A continuous wavelet converts a time-series signal into the time–frequency domain. It depicts the magnitude and periodicity of the time-series data [X_t] in the time domain (Mahasa et al., 2023). The Morlet wavelet transform of a time series (X_t) is defined as the convolution of the series with a set of “wavelet daughters” generated by the mother wavelet by time translation and scale factors (8): W a v e τ , s = ∑ t x t 1 √ s Ψ * t − τ s , (8)

where s is the wavelet scale, τ defines the position of the wavelet window in the time or the translated time index, and Ψ is the mother Morlet function with * denoting the complex conjugate. In the present study, the Morlet wavelet function (Ψ) as a basis of the wavelet function is used through the following expression (9): Ψ t = π − 1 / 4 e i ω t e − t 2 / 2 , (9)

where ω is the rotation rate or angular frequency per unit time in radians and t is the time step. One revolution is equal to 2π (radians); therefore, the period (or inverse frequency) measured in time units equals 2π/6. The wavelet power can be expressed through the following expression (10): P o w e r τ , s = 1 s W a v e τ , s 2 . (10)

A cross-wavelet transformation (CWT) function was applied to measure the relationship between NDVI and individual climatic variables using a similar time step. The cross-wavelet [W_xy (τ,s)] is expressed in the following way (11): W x y τ , s = w x   τ , s . w y *   τ , s , (11)

where * denotes the complex conjugate of the y time series. The following expressions can be used for the determination of the power of cross-wavelet, which is the measure of coherence between pairs of time-series data (12 and 13): C o h e r e n c e = s . W a v e . x y 2 s P o w e r . x . s P o w e r . y , (12) P o w e r . x = 1 s W a v e τ , s 2 ; P o w e r . y = 1 s W a v e τ , s 2 . (13)

The difference in phase between two time-series data depicts an in-phase oscillation pattern of two time-series data if they are matched with each other; in contrast, the difference in phase depicts an out-of-phase pattern if the two time-series data do not match (Grinsted et al., 2004). The lagging or the leading tendencies of the variables are indicated through the arrows (Schmidbauer and Roesch, 2018). The coherence between the variables was evaluated at the p-value of 0.05 and was incorporated in the cone of influence (Pal and Devara, 2012). The entire process of wavelet analysis and its graphical expression was performed in RStudio using the WaveletComp package (Schmidbauer and Roesch, 2018). To detect the presence or absence of multicollinearity among climatic variables, the variance inflation factor (VIF) was obtained in RStudio using the faraway package (Faraway, 2022). The spatial distribution of NDVI was depicted in the Google Earth Engine (GEE) code editor. The Landsat 5 TM dataset was used for the depiction of the year 2000 NDVI data, and the Landsat 7 ETM+ dataset was adopted for the recent year. In the cases of 2000 and 2021, the mean annual time frame was considered for the depiction of NDVI. Finally, the results were downloaded and represented in ArcGIS 10.5.

The relationship of mean NDVI to the climatic factors was determined by using Karl Pearson’s coefficient of correlation. Precipitation showed a positive and moderate correlation with mean NDVI (r = +0.50). Soil moisture and precipitation are important climatic factors that have a direct impact on the growth of plants and crops (Chen et al., 2023; Yang et al., 2023). Evapotranspiration depicted a higher positive correlation with all the factors including mean NDVI (r = +0.71). On the other hand, a relatively insignificant correlation was observed between the 2-m air temperature and mean NDVI (r = +0.03). Soil moisture and NDVI depicted a high positive correlation of +0.81 (Figure 2). The positive correlation between the mean NDVI and climatic factors depicted the definite influence of climatic factors on the growth of plants and crops (Tian et al., 2019; Wang et al., 2022).

Before further assessing the relationship between climatic factors with NDVI, it is crucial to evaluate if any kind of multicollinearity exists between climatic factors (Li M et al., 2021; Zhang et al., 2021). The presence of multicollinearity can result in unstable estimations as it enhances the variance of the coefficient of regression (Chan et al., 2022; Zhu et al., 2023). VIF is a widely used method adopted by scholars of different disciplines for depicting multicollinearity among variables (Chan et al., 2022; Park et al., 2022). The VIFs for the climatic factors were 2.43, 4.85, 7.06, and 2.12 for precipitation, soil moisture, evapotranspiration, and 2-m air temperature, respectively, with a p-value of <2.2×10⁻¹⁶. All the VIFs showed values less than 10 with a very low p-value, which indicated the absence of any significant collinearity among climatic factors (Akinwande et al., 2015).

Figure 3 depicts the spatial dynamics of NDVI from 2000 to 2021. The illustration clearly shows a significant change in vegetation cover over the last 2 decades. In recent years, the greener segments, which are mostly associated with agricultural activities, have expanded in the northern and eastern segments. In previous years, the concentration of crops/vegetation was mostly observed in the central segment.

The linear trend of mean NDVI showed a clear upward trend (Figure 4A) for the period 2000–2021. The range of mean NDVI from 2000–2021 was between 0.12 and 0.49 as minimum and maximum, respectively, with a mean of 0.07. The monthly minimum NDVI value was observed during May 2000 (pre-monsoon), while the maximum NDVI value was observed during September 2011 (monsoon). The standard deviation of 0.08 showed relatively a small deviation in the entire time span with a positive skewness (+1.17). MK trend analysis depicted a positive trend at the 95% confidence level with a Sen’s slope of 3.3×10⁻⁴.

The precipitation rate varies between 0 mm/day and 287.41 mm/day in the entire time span. The maximum precipitation was observed during August 2012 (monsoon). The mean precipitation value was as high as 46.17 mm/day, the standard deviation was as high as 63.18 mm/day, and skewness indicated a high positive value (+2.06). The MK trend test and Sen’s slope indicated a downward trend in the time frame (z = −2.46 and Sen’s slope −6.88×10⁻³, respectively) (Figure 4B). Soil moisture ranged between 47.95 kg m^–2 and 86.97 kg m^–2 with a mean of 6.52 kg m^-–2 (Figure 4C). The standard deviation was as high as 8.29 kg m^–2, while skewness was 0.71. The MK trend test and Sen’s slope indicated an upward trend in the time frame (z = 2.71 and Sen’s slope 1.72×10⁻², respectively). Evapotranspiration ranged between 8.73×10⁻⁷ kg m^–2 and 4.22×10⁻⁵ kg m^–2 with a mean of 8.88×10⁻⁶ kg m^–2 (Figure 4D). The standard deviation and skewness were 1.09×10⁻⁵ kg m^–2 and 1.09, respectively, while the MK trend and Sen’s slope were 1.42×10⁻⁸. Figure 4E shows the trend of 2-m air temperature on NDVI. The range of 2-m air temperature from 2000 to 2021 was 12.52°C (minimum) and 37.21°C (maximum) (Figure 4E). The standard deviation was high, at 7.30, while the skewness was −0.34. The minimum temperature was observed during January 2020 (winter), and the maximum temperature was observed during June 2002 (pre-monsoon). MK trend statistics showed a downward trend (z = −1.16) with a Sen’s slope of−5.66×10⁻³.

The mean monthly characteristics of NDVI depicted some important facts and complex relationships with other climatic factors. The months of May and June showed the lowest means (0.14 and 0.15, respectively) and minimum inter-quartile ranges (IQRs) (0.02 and 0.02, respectively) (Figure 5A). On the other hand, August and September depicted the highest means (0.37 and 0.36, respectively) and highest IQRs (0.13 and 0.13, respectively). The maximum NDVI values can be observed in July, August, and September. Precipitation (Figure 5B), soil moisture (Figure 5C), and evapotranspiration (Figure 5D) confirm a similar pattern; 2-m air temperature is an exception (Figure 5E). During June, 2-m air temperature reached its zenith with a mean and IQR of 35.54 and 1.75, respectively.

Figure 5E depicts the mean monthly temperature pattern of the entire time span. In the selected study area, the winter season spans between mid-November and February, and the pre-monsoon season or summer prevails from April to June. The monsoon season ranges between July and September. The characteristics of the post-monsoon season can only be observed during October and early November (Dhakar et al., 2013). Seasonal characteristics were depicted by observing the behavior of all the variables in the four major seasons, viz., winter, pre-monsoon, monsoon, and post-monsoon. During winter and pre-monsoon, NDVI showed relatively low IQRs (0.04 and 0.05, respectively) (Figure 6A) (Table 2) with a very low standard deviation (0.03) and depicted a significant association with soil moisture in particular (Figure 6B). During the monsoon season, IQR showed the highest value of 0.20 and depicted the maximum association with the climatic factors. During monsoon, the IQRs of related climatic factors like precipitation (Figure 6C), soil moisture (Figure 6D), and evapotranspiration (Figure 6C) were 119.69 mm/day, 16.41 kgm^–2, and 1.7×10⁻⁵ kgm^–2, respectively. At the same time, the IQR of the 2-m air temperature was at a minimum during the monsoon season (3.51) (Figure 6E). In post-monsoon, the IQR of the mean NDVI was 0.13, while at the same time, the IQR of precipitation was 44.80 with a considerably lower mean value. Soil moisture also follows a similar pattern as precipitation in terms of IQR. In the case of NDVI, the standard deviation was relatively smaller in all seasons except monsoon. In contrast, all of the climatic variables showed a higher standard deviation. Table 2 summarizes the descriptive statistics of NDVI and all the climatic factors.

The precipitation time series and NDVI time series depicted some significant results. A continuous segment of coherency among the factors can be observed in the eighth to the 16th month period. The highest cross-wavelet power (CWP) level can be observed after 2010 in the 12th month segment period (Figure 7A). Another discrete segment of coherency was observed in the fourth to the seventh month period. The direction of the arrows confirms that in these segments, the factors were associated with an in-phase condition. Results obtained from the cross-wavelet analysis of evapotranspiration and NDVI depicted some significant facts. The CWP level ranged between below 0.3 to above 1.5. A continuous segment of coherency with a longer periodicity of the evapotranspiration time series and NDVI time series can be observed from the ninth to the 12th month period throughout the temporal dimension. A segment of periodicity can also be observed between the fifth and the seventh month period (Figure 7B). In this case, the segment was discrete before 2010, while it is relatively continuous afterward. The maximum CWP level was observed after 2010 and in the period segment of the 11th to the 12th month period in a discrete manner. The arrows in the segments toward the right (fifth to the seventh month period) and top-right (ninth to the 12th month period) directions indicate the in-phase cyclicity of evapotranspiration and NDVI (Schmidbauer and Roesch, 2018; Liu et al., 2023). In the case of the ninth to the 12th month period, the segment was associated with the situation when the mean NDVI was leading and evapotranspiration was lagging. Soil moisture and NDVI showed a similar pattern of coherency and phase condition except for CWP level. In this context, most of the area was associated with a relatively lower CWP level, with only a few patches with a higher CWP level. An in-phase condition can be observed during the fourth to seventh month period and the ninth to the 16th month period (Figure 7C). In the case of 2-m air temperate and NDVI, only a single segment of coherency can be observed in the eighth to the 16th month period with the highest CWP level between 2010 and 2015. The direction of the arrow indicates a slight in-phase condition (Figure 7D).