## Abstract

Beside *in situ* observations, satellite-based products can provide an ideal data source for spatiotemporal monitoring of drought. In this study, the spatiotemporal pattern of drought was investigated for the northwest part of Iran using ground- and satellite-based datasets. First, the Standardized Precipitation Index series were calculated via precipitation data of 29 sites located in the selected area and the CPC Merged Analysis of Precipitation satellite. The Maximal Overlap Discrete Wavelet Transform (MODWT) was used for obtaining the temporal features of time series, and further decomposition was performed using Ensemble Empirical Mode Decomposition (EEMD) to have more stationary time series. Then, multiscale zoning was done based on subseries energy values via two clustering methods, namely the self-organizing map and K-means. The results showed that the MODWT–EEMD–K-means method successfully identified homogenous drought areas. On the other hand, correlation between the satellite sensor data (i.e. the Normalized Difference Vegetation Index, the Vegetation Condition Index, the Vegetation Healthy Index, and the Temperature Condition Index) was evaluated. The possible links between central stations of clusters and satellite-based indices were assessed via the wavelet coherence method. The results revealed that all applied satellite-based indices had significant statistical correlations with the ground-based drought index within a certain period.

## HIGHLIGHTS

Discussing the spatiotemporal variations of drought using

*in situ*observations and satellite-based datasets.Applying the hybrid MODWT–EEMD–K-means method for catching similar zones.

Discussing the possible links between the SPIs of the central stations of clusters and satellite-based drought indices via the wavelet coherence (WTC) method.

Evaluating drought conditions using the satellite-derived NDVI and LST products.

### Graphical Abstract

## INTRODUCTION

Drought is considered as a natural disaster, effects of which accumulate slowly over a long period. During the past decades, frequent and severe droughts caused considerable ecological, agricultural, and financial loss around the world (Van Loon 2015; Saghafian & Sanginabadi 2020). Drought is generally difficult to define and understand due to its invidious nature (Vicente-Serrano 2007; Dadson *et al.* 2019; Li *et al.* 2020). The type and severity of drought vary spatially, and the resulting impacts depend on the vulnerability of affected sectors (Zhou *et al.* 2018; Kisi *et al.* 2019). Monitoring the spatiotemporal variation of drought makes it possible to reduce vulnerability and establish adaptation strategies (Nischitha *et al.* 2014). There are different techniques for quantitative analysis of drought phenomenon, and among them drought indices have been widely applied for this aim. For instance, Zhou & Zhou (2021) investigated drought duration, severity, intensity, affected area, and centroids during 1960–2015 for the Pearl River basin based on the Standardized Evapotranspiration Deficit Index and three-dimensional clustering algorithm and then revealed how these drought characteristics have affected net primary productivity. Li *et al.* (2020) used the Standardized Precipitation Evaporation Index time series for classifying drought grades in the Guanzhong area based on the meteorological and remote sensing datasets. Among drought indices, the Standardized Precipitation Index (SPI) is an ideal tool in drought severity characterizing (Li *et al.* 2015). Nowadays, this index is known as one of the most common methods for drought monitoring (Kwak *et al.* 2016). The SPI index has statistical consistency and is capable of showing short- and long-term drought impacts on precipitation anomalies. Although the drought indices provide useful information in drought assessment, they are derived from gauge observations and have several limitations, such as low spatial resolution and discontinuity. The satellite-based products can provide alternative datasets with higher spatial resolution, which are temporally and spatially continuous (Yilmaz *et al.* 2008). Therefore, monitoring of drought conditions is possible in ungauged or inaccessible areas via the satellite datasets (Skakun *et al.* 2016). Although the satellite-based products represented have been widely used in hydrological simulations, raw satellite datasets have some disadvantages, such as the amount of cloudiness which affects the accuracy of the results. Also, the data quality of these products will inevitably be affected by observational bias, spatial scale, and retrieval method. Therefore, satellite-based datasets need to be examined and verified by comparison with ground-based data. Satellite products have been shown to play an important role in assessing drought-related vegetation conditions, especially in semi-arid and arid regions. Owing to the close relationship between available soil moisture and vegetation vigor, two indices, the Normalized Difference Vegetation Index (NDVI) and the Land Surface Temperature (LST) index, have been widely used to assess drought conditions. The NDVI is an indicator of green biomass, leaf area index, and pattern of production. According to Zhou *et al.* (2012), the NDVI principle is that the internal mesophyll structure of healthy green leaves strongly reflects near infrared radiation and leaf chlorophyll, and other pigments absorb a large proportion of the red visible red radiation. This becomes reverse in the case of unhealthy or water-stressed vegetation. Also, the LST is used to show the energy balance at the surface of earth, one of the important variables in the physics of land-surface processes on regional and global scales.

On the contrary, within the past few decades, the time-series pre-processing approaches are commonly applied for hydrometeorological data analyses, given the seasonal fluctuations and non-stationary features of signals (Roushangar *et al.* 2021a). The Discrete Wavelet Transform (DWT) is a common approach used by many researchers for time-series decomposition (Agarwal *et al.* 2016). The modified version of the DWT is termed as the Maximal Overlap Discrete Wavelet Transform (MODWT). In this latter approach, time series should be projected onto a set of non-orthogonal basis functions in order to generate the wavelet coefficients. Another technique used for signal decomposition is the Ensemble Empirical Mode Decomposition (EEMD) method. According to Rezaie-Balf *et al.* (2019), this method is an intuitive, empirical, and self-adaptive data processing tool. The EEMD method is especially appropriate for analysis of nonlinear and non-stationary time series, which comes out as a superiority of the EEMD. In addition, clustering methods are used to identify homogenous drought regions. Clustering can be defined as the classification of a large number of data into a smaller number of groups so that data belonging to the same cluster are as similar as possible, while they have heterogeneity with the data belonging to other groups (Barton *et al.* 2016). The clustering approaches regard the spatiotemporal domains of components; so, the obtained outcomes from them indicate the local characteristics and temporary variations. The K-means and self-organizing map (SOM) methods are among the most common clustering methods, which were adopted in the present study.

Drought is a quickly increasing risk that can lead to hardly repairable or unrepairable damages to the nature as well as socio-economy, especially for arid and semi-arid regions. Analyzing the spatiotemporal pattern of drought can provide a better understanding of drought characteristics in these regions. Meteorological data mainly originate from a ground station network, and the quality and stability of them cannot be guaranteed due to the scarcity of stations, uneven spatial distributions, limited time scales, and vulnerability to environmental and human factors. Since meteorological data, such as precipitation and land air surface temperature, collected by surface observation stations often possess poor spatial resolution, especially in remote regions with difficult access and in some developing countries, satellite-based data and the remotely retrieved NDVI and LST data may provide a valuable source of information for monitoring drought. Therefore, in this study, a multi-spatiotemporal method based on both *in situ* observations and satellite datasets was proposed to analyze the drought correlations with the satellite sensor datasets in the northwest Iran. In this regard, first, the 3-month SPI series were calculated via (1) precipitation data obtained from 29 sites located in the selected area and (2) the CPC Merged Analysis of Precipitation (CMAP) satellite covering the period of 1990–2020. In the second step, the MODWT decomposition was applied for obtaining the temporal features of the signals, and further decomposition was performed via the EEMD. This step is generally termed as the pre-processing of the data. The time-series decomposition into various periodicity scales with the serial application of MODWT and EEMD results in more stationary subseries. Then, the energy values of each subseries were computed, and the mean energy values of each station were used for spatial categorization (i.e. clustering) of the area. Two clustering techniques, namely the K-means and SOM, were used for this aim.

In the next step, three satellite-derived NDVI and LST products, namely the Vegetation Condition Index (VCI), the Temperature Condition Index (TCI), and the Vegetation Healthy Index (VHI), were used to evaluate drought conditions. The NDVI and LST data were obtained from the MODIS satellite datasets MOD13Q1 and MOD11A2, respectively. The Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data series used for the present study covered the period of 2000–2020. From the remote sensing technologies, MODIS, Long-term Data Record, and Advanced Very High Resolution Radiometer (AVHRR) datasets have been extensively used to assess vegetation conditions and drought all over the world (El-Vilaly *et al.* 2018). Among these datasets, MODIS is more popular because of its heritage sensors in terms of spatial, spectral, and temporal resolutions. A major advantage of using the MODIS imagery is the availability of a suite of products ranging from raw images to highly processed products such as vegetation indices. The advantage of MODIS data over National Oceanic and Atmospheric Administration (NOAA)–AVHRR data is that the MODIS data provide continuity for historical applications of time series. The MODIS NDVI data are computed from atmospherically corrected bi-directional surface reflectance that have been masked for water, clouds, heavy aerosols, and cloud shadow. Therefore, in this study, the MODIS satellite datasets were used for calculating the NDVI and LST time series. On the other hand, the appropriate correlation of precipitation data of the CMAP satellite with the ground-based data has been proved in a previous study of the authors (Roushangar *et al.* 2021b). Finally, for further evaluating the relationships of SPIs with the NDVI and LST products, any possible links between the SPIs of the central stations of clusters and satellite-based indices were investigated using the wavelet coherence (WTC) method. To the best of the authors’ knowledge, there are no known studies dealing with the relationships between the SPIs and satellite-based indices in the selected area. The novelty of this study is that for the first time, we are able to determine the statistical relationships and possible links among the applied drought indices in the northwest Iran. This study can improve the insight into the relationships between the SPI and the NDVI, LST, VCI, and TCI, which helps decision-makers manage droughts and provide guidance for drought-related studies across other regions.

## MATERIALS AND METHODS

### Proposed methodology

In this study, the spatiotemporal variations of drought were assessed for the northwest Iran using both *in situ* observation and remote sensing-based datasets. As shown in Figure 1, the procedures of the multi-spatiotemporal model could be summarized as follows:

Obtaining the precipitation data from the

*in situ*stations and the CMAP satellite during the period of 1990–2020.Calculating the SPIs time series.

De-noising the SPIs time series.

Decomposition of the de-noised series into several subseries via the MODWT and further decomposition of subseries with the EEMD (the pre-processing).

Calculating the energy values of each subseries.

Zoning the selected region based on the energy amounts of the subseries and the mean SPI values and investigating the relationship between them.

Spatially identifying the homogenous drought regions via the K-means and SOM clustering techniques.

Validating the outcome of spatial clustering methods using the following three indices: the Davies–Bouldin Index (DBi), the Dunn Index (DI), and the Silhouette coefficient (SC).

Selecting the central station of each cluster.

Obtaining the NDVI and LST time series from the MODIS satellite covering the period of 2000–2020.

Calculating the VCI, TCI, and VHI indices from the NDVI and LST time series.

Determining the possible links between the central stations of clusters (ground-based data) and satellite-based drought indices. The applied methodology is explained in more detail in the following sections.

### Description of the study area

Owing to the high variability of annual precipitation in the northwest part of Iran, drought is more frequent in this area; therefore, this area is selected for assessing the spatiotemporal variation of drought. The precipitation datasets of 29 sites obtained from the meteorology portal of Iran (http://www.weather.ir) were used for this aim. Beside the *in situ* observation data, the CMAP remote sensing product was employed to compute the SPI time series for the selected stations. According to Xie & Arkin (1996), the CMAP product is a collection of precipitation dataset, which is generated using remote sensing-based precipitation estimations and gauge data analysis. Also, the NDVI and LST series were obtained from the MODIS on the National Aeronautics and Space Administration (NASA) Earth Observation Satellite (EOS) Terra (https://modis.gsfc.nasa.gov/). Figure 2 shows the geographical locations of selected sites.

### Pre-processing approaches

*et al.*2006; Roushangar

*et al.*2018). This process enhances the quality of the time series. The MODWT is one of these techniques, which is essentially the modified version of the DWT as stated above. The MODWT can be applied for decomposition of original time series into an approximation component (scaling coefficient) and many detail components (wavelet coefficients). The scaling coefficients contain the low-frequency information that is the most crucial section in providing the signal identity, whereas the wavelet coefficients reveal flavor of the signal (Sujjaviriyasup 2017). The expression for the MODWT is as follows:where

*X*is the time series,

*W*is the wavelet details,

_{j}*T*is the overall trend of original signal, and

_{N}*L*is the number of decomposition level. According to Dghais & Ismail (2013), the MODWT is non-orthonormal and shift-invariant, and it is applied for time series with different sizes. For more details on the MODWT, the reader is referred to Percival & Walden (2000).

The empirical mode decomposition (EMD) is another efficient method for pre-processing of signals. In the EMD, the signal is broken down into several inherent mode functions (IMFs), which makes it possible to process nonlinear time series. The capability to specify the instantaneous frequency of the signal is one of the benefits of EMD. The EEMD is the modified version of the EMD. According to Wu & Huang (2009), the EEMD is used for solving the EMD mode mixing issue. Both the EEMD and the MODWT facilitate incredible vision in frequency and time domain of data to capture the seasonal and nonlinear properties.

*E*) of the subseries in order to reduce the number of inputs in the area clustering process.where SP shows the amount of SPIs decomposed via the MODWT and EEMD, and

*m*is the SPIs month.

### Area clustering techniques

Clustering is the task of dividing a set of patterns or data points into a number of groups, such that data points in the same groups are more similar to other data points in the same group than those in other groups. In simple words, the aim is to segregate groups with similar traits and assign them into clusters. Among different clustering techniques, the K-means generally applied for data-mining cluster analysis and vector quantization (Pham *et al.* 2005). Using the K-means, *n* observations are classified into *k* clusters and in this process, each entity (observation) is assigned to the cluster with the nearest mean (Pham *et al.* 2005). Based on Likas *et al.* (2003), the following steps should be performed in data clustering through the K-means method:

Selecting an initial number of clusters and defining the centers;

Assigning each observation to a cluster with the closest center considering its Euclidean distance;

Reordering the positions of the centers after assigning all observation to one cluster;

Repeating steps 2 and 3 until the memberships of cluster do not vary.

Given the ability of K-means in clustering large datasets as well as its implementation in a short time, it is one of the most popular clustering methods. Additionally, the SOM is another commonly used approach in exploring and extracting the inter-relationships of higher dimensional datasets. This method is useful in data clustering and anticipating in a widespread range of fields (Cereghino & Park 2009). The SOM has high ability in extracting implicit patterns from high-dimensional dataset and grouping the captured patterns into a low-dimensional output layer (Hsu & Li 2010). In this process, similar inputs remain close together in the output neurons, while the structure of data is preserved. Since the neurons of the output layer are usually ordered in two-dimensional grids, the produced topology can be visualized to provide an insight into the considered system (Kalteh *et al.* 2008; Iwashita *et al.* 2018).

The outcomes of spatial clustering were evaluated based on the following three validation metrics: DI, DBi, and SC. The DI goal is to distinguish a category of clusters that are well-set, with a small variance among components of the cluster, and well-detached, where the averages of the various clusters are adequately far apart when compared to the within-cluster variance. A higher DI value shows a better clustering outcome as it shows a well-compacted cluster (Agarwal *et al.* 2016). Computational cost increases when the number of clusters and dimensionality increases, which is a disadvantage for the DI. The DBi is a widely applied internal evaluation criterion, which is applied to distinguish the number of optimal clusters that are well-detached and well-set based on the content and specification of dataset. A lower DBi value represents better clustering results. On the other hand, the DBi has a disadvantage in that best information detection cannot be implied by a good reported DBi value (Kasturi *et al.* 2003). Also, the values of SC indicate the degree of similarity for samples within the cluster. The range of SC varies between −1 and 1, and higher values indicate that the samples are well-matched to their own clusters (Hsu & Li 2010).

### WTC approach

*S*is a smoothing operator and is the cross-wavelet transform. This deﬁnition is very similar to that of the traditional correlation coefficient. It is beneficial to consider the WTC as a localized correlation coefficient in a time–frequency space (Chang

*et al.*2019). The

*S*operator can be written as:where

*S*

_{scale}and

*S*

_{time}signify smoothing along the wavelet scale axis and in time, respectively. An appropriate smoothing operator for the Morlet wavelet can be expressed as given in Equations (5) and (6), in which

*Π*is the rectangle function,

*c*

_{1}and

*c*

_{2}are normalization constants, and the factor of 0.6 is the empirically determined scale de-correlation length for the Morlet wavelet (Grinsted

*et al.*2004).

*c*

_{1}and

*c*

_{2}are specified numerically after the calculation of the two convolutions (the ﬁlters

*c*

_{1}Π (0.6

*S*) and must have unit weight):

## RESULTS AND DISCUSSION

### Decomposition of the SPIs obtained from the *in situ* observations and remote sensing data

For assessing the spatiotemporal variations of the drought time series in the considered region, first, the MODWT and EEMD methods were applied to decompose the de-noised SPIs into several subseries. Since different internal and external factors may affect the quality of data during the data collection process, the wavelet transform was applied to de-noise the original SPIs datasets. To decompose a time series via the MODWT, at first, an appropriate wavelet function (a function with similarities to the original signal as much as possible) should be selected. Therefore, different wavelet functions (i.e. Daubechies, Symlets, Haar, Coiflets, Meyer, Bior, and Rbio) were tested. On the other hand, the optimum decomposition level (*L*) was determined via a trial-and-error process. For this purpose, the best wavelet filters for *L* = 4 were determined and then, by changing the value of *L* parameter, its optimized amount was sought. The root-mean-square error criterion was used to select the best case. From Figure 3, it is seen that the db4 resulted in the best outcome among the applied wavelet filters. Also, the optimized decomposition level was found to be *L* = 4. Khan *et al.* (2020) used wavelet-based hybrid models for meteorological drought modeling and stated that among different mother wavelets, the debauches (*db*) type led to better results. This wavelet function has full scaling and translational orthonormality properties and also signifies that the wavelets have non-zero basis functions over a finite interval.

After signal decomposition, five subseries (i.e. four wavelet coefficients (*W _{i}*) and one scaling coefficient (

*V*

_{4})) were extracted. In Figure 4, the decomposed subseries of the SPIs obtained from the

*in situ*observations and CMAP satellite data are shown for the S18 station. As seen in this figure, the lower wavelet levels expressed higher frequency components, while the higher levels showed components with lower frequencies (including the trend of signal). Each wavelet coefficients represent a specific time period. For example, for monthly datasets, the

*W*

_{1},

*W*

_{2},

*W*

_{3}, and

*W*

_{4}show 2

^{1}= 2, 2

^{2}= 4, 2

^{3}= 8, and 2

^{4}= 16-month periods, respectively. In the second step, the first two subseries of wavelet coefficients (i.e.

*W*

_{1}and

*W*

_{2}) were further decomposed via the EEMD, in order to have more stationary subseries. Using the EEMD, the signal is decomposed into several IMFs and one residual subseries. The sum of these subseries will be the same original signal. Owing to the increase in the number of inputs after signal decomposition, for decreasing the number of inputs, the energy values of subseries were computed and applied as inputs of the clustering techniques.

In the following step, the selected region was zoned via the obtained mean energy values to investigate the relationship between energy values and the SPIs time series. Figure 5 shows the results of the area zoning using the mean energy values and the mean drought index for both *in situ* observations and the CMAP satellite datasets. It should be noted that in the zoning process, the inverse distance weight (IDW) was used for spatial interpolation. The interpolation methods can be classified into the following two major groups: deterministic and geostatistical. Deterministic interpolation techniques create surfaces from measured points, based on either the extent of similarity or the degree of smoothing. A deterministic interpolation can either force the resulting surface to pass through the data values or not. On the other hand, geostatistical (or stochastic) methods capitalize on the spatial correlation between neighboring observations to predict attributed values at unsampled locations. The IDW is a determinist method and estimates the values of an attribute at unsampled points using a linear combination of values at sampled points weighted by an inverse function of the distance from the point of interest to the sampled points. The assumption is that sampled points closer to the unsampled points are more similar to it than those located further away. This method has been widely used in many different fields, such as hydrology and earth science (Jeong *et al.* 2020). According to Figure 5, for both datasets, the southeast part suffered from severe drought, while drought was slighter in the southwest and northwest parts of the study area. In the southeast part of the selected region, the SPI amounts were higher than the other parts, and in the southwest part the SPI values reduced and reached the lowest value. The lowest energy values were obtained for the central parts, where there were more arid regions. Compact counter lines were observed on the northwest parts of the selected area, meaning that energy amounts varied rapidly on these parts which were located on rainy regions of the area. Spatial changes of energy became smoother for the eastern zones which were semi-arid regions. From these results, it is clear that there is an inverse relationship between the drought index and energy values, i.e. the smaller the energy values, the higher the intensity of the drought index is.

### Results of the SPI regionalization applied to the *in situ* observations and satellite data

To cluster 29 stations in the selected area, the mean energy values of each station were used as inputs of the K-means clustering technique. According to Figure 6, three SC, DBI, and DI evaluation criteria were applied for determining the optimal number of clusters (NC). The DI is defined as the ratio of the minimal intra-cluster distance to the maximal inter-cluster distance. This index is an internal evaluation scheme for evaluating clustering algorithms, where the result is based on the clustered data itself. A large value of the DI is preferred as it represents a better separation. The DBi determines an optimal cluster number, such that clusters are well separated and compact. Since the objective is to obtain clusters with minimum intra-cluster distances, a small DBi value is important. Also, if most members have a high SC value, then the formation of the clustering is suitable. On the contrary, if many members have a low or negative SC value, accordingly the clustering formation may have too many or too few clusters. The number of clusters to adequately represent the dynamic properties of drought data was determined to be between 2 and 8. Following the application of the K-means clustering technique with a training program which ran through 1,000 trials, the results obtained for various numbers of energy-based clusters were compared and the optimal number of clusters (Figure 6) was determined through the comparison of three validity indices. In the clustering process via the K-means approach, it could be seen that for the *in situ* dataset, a cluster number of 4 (–) yielded values of 0.89, 1.36, and 1.82 for the SC, DBI, and DI, respectively, and for the CMAP data, the NC = 4 yielded values of 0.77, 1.36, and 1.76 for the SC, DBI, and DI, respectively, indicating a better performance in determining the homogenous regions compared with other NC values. Therefore, the NC = 4 was selected as being the optimum value for the clustering of stations. Also, the SOM method was applied for clustering the 29 stations into a visible two-dimensional topology map. In this regard, different map sizes were tested (i.e. 2 × 2 to 10 × 10). The results indicated that in the case of *in situ* dataset, the NC = 5 led to better outcomes in which values for the SC, DBI, and DI were 0.78, 1.44, and 1.78, respectively. While, in the case of satellite dataset, the NC = 6 performed better in determining homogenous regions due to having 0.68, 1.45, and 1.57 values for the SC, DBI, and DI indices, respectively. However, based on the results, a failure has occurred in the SOM outcome with NC = 5. In this case, there were three clusters with only one station as membership of the cluster, and two clusters with >10 memberships, which is an undesired result. Considering both performance criteria and classification of stations outcomes, it was proved that the K-means performed better than the SOM; thus, the K-means method was selected for further analysis.

Decomposed subseries in the pre-processing step show different monthly scales; however, some of these components have no strong correlation with the original SPIs. Therefore, mean energy values of the subseries were computed and considered as inputs of the K-means for performing spatial categorizing. The dynamic properties of time series based on their energy values can be used to improve the efficiency of the clustering approaches. The geographical location of stations based on clustering through energy values is shown in Figure 7(a). In Figure 7(b), the members of each cluster and their center stations are presented based on the SC index. From Figure 7, it can be seen that some of the stations located in the same cluster are spread across the study area, indicating that the geographic contiguity is not the sole basis of clustering given the differences in energy values calculated for each station.

### Correlation analysis between drought indices obtained from *in situ* observations and remote sensing datasets

*et al.*(2018), the NDVI not only shows the presence of vegetation on a pixel basis, but also provides measures of the condition or amount of vegetation within a pixel. Also, the LST is regarded as a good indicator of the energy balance at the surface of Earth due to its importance in the physics of land-surface processes on both regional and global scales (Zhang

*et al.*2015). In this study, the MODIS (MOD13Q1) and MODIS (MOD11A2) datasets were used for obtaining the NDVI and LST time series, respectively. In addition, the study investigated the potential use of the TCI and VCI data extracted from the multi-temporal LST and the NDVI. According to Zhang

*et al.*(2015), drought monitoring through only vegetation condition indicators has some limitations. However, they can be solved via incorporating surface temperature (LST). So, the VHI can be used as an appropriate alternative in monitoring of drought. The mentioned indices were calculated as:where NDVI

*and LST*

_{i}*show the NDVI and LST values of*

_{i}*i*th month, respectively. The NDVI

_{max}and NDVI

_{min}are the multi-year (2000–2020) absolute maximum and minimum NDVI amounts. The LST

_{min}and LST

_{max}show the smoothed multi-year minimum and maximum LST for every pixel at a specific period. Also,

*α*= 0.5 indicates equal contribution of the TCI and VCI. The obtained results are shown in Figures 8 and 9 and are listed in Table 1. Figure 8 shows the spatial distribution of the NDVI for a dry (2008) and a wet (2010) year, and also for the mean NDVI during the considered period. According to this figure, the southwest part of the study area was less affected by drought, while severe droughts occurred in the western regions. From the results, the long-term (2000–2020) mean NDVI values in the selected area were between 0.008 and 0.45, while the LST values were between 26.2 and 31.68 °C. In general, a region with a lower NDVI value indicates a high LST value and vice versa. Based on Guha

*et al.*(2017), the NDVI with negative values indicate water bodies, the NDVI with positive values up to 0.2 designates bare areas, and the NDVI values >0.2 indicates vegetation. The higher NDVI value shows the healthier vegetation with high chlorophyll content. Also, a region with a low VHI indicates the high occurrence of drought, while a higher VHI value shows wet conditions. From the VHI time series, it can be inferred that the moderate (20 < VHI < 30) and mild (30 < VHI < 40) droughts occurred in the region. The results of the relationship between the LST and the NDVI showed that there was a negative correlation between the two mentioned indices. By increasing the LST to a certain amount, the NDVI reduced since the LST was stronger during the day than at night. Zhang & Jia (2013) also reported similar finding. They reported that the LST is stronger during the daytime and this causes a higher level of vegetation stress that can lead to the incidence of drought. Also, Karnieli

*et al.*(2010) investigated the use of the NDVI and LST for drought assessment and indicated that the increases in evaporation along with a decrease in soil moisture caused by higher temperatures affect the NDVI. Based on Gidey

*et al.*(2018), during drought, only the LST increased while the NDVI significantly decreased due to poor moisture availability and rising surface temperature. This can be a result of low precipitation in their study area. However, in low-temperature areas, the relationship between the NDVI and the LST might be positive. For instance, during the winter or cold season, the regression between the NDVI and LST can be positive (Sun & Kafatos 2007). A negative relationship was also observed between the TCI and the VCI. The results showed a positive relationship between the VHI and the SPI. Gidey

*et al.*(2018) indicated that when the SPI increases, the VHI also tends to increase. The statistical relationship between the SPI and the VHI depicted how the drought monitored by the VHI associates with the meteorological drought measured by the SPI. The regression analysis result showed that the relationship between the SPI and the VHI in the 3-month timescale was positive, strong and statistically significant. Also, Yan

*et al.*(2016) stated that the statistical relationship between the SPI and the VHI has strong correlation during a drought period in southwest China.

### Analysis of the relationship between the center stations of each cluster with the satellite drought indices

In this section, the WTC analysis was applied for assessing the relationships between *in situ* observations and satellite-based drought indices in the study area. The WTC was used to find out whether there was a significant oscillatory period between the selected drought indices. Since there were 29 stations in the selected region, the use of wavelet dependencies on all stations and all variables was not possible; so, the central stations of each cluster were considered as representative of all stations. After clustering the area, a number of stations are assigned to each cluster as memberships of that cluster, and one of them is the central station. The central station should have the lowest distance from the cluster's center and also have the highest SC value compared to other stations. These members are marked in yellow in Figure 7(b) for all clusters. Finally, four stations, namely S12, S13, S10, and S21, were selected as representatives of each cluster. The results of WTC analysis between the drought indices are shown in Figure 10. A relative phasing of two time series is provided on each of the panel in this figure via arrows for indicating the cause and effect relations between applied indices. Cold parts outside the signiﬁcant areas (bluish regions) showed frequencies and time with no dependence in the time series. The causality relation between two time series can be assessed using the phase patterns. The direction of the importance of arrows can be described as the following:

*X*(*t*) (SPI) and*Z*(*t*) (second time series, i.e. VCI, TCI, or VHI) are in-phase (positively related) when arrows point to the right.*X*(*t*) and*Z*(*t*) will be anti-phase (negatively related) when the direction of the arrows is to the left.Arrows pointing in other directions show leads or lags between two time series. For instance, arrows with the straight down direction denote that

*Z*(*t*) leads*X*(*t*), while arrows pointing up signify that*X*(*t*) leads*Z*(*t*).

For having a correct sight about the WTC graph, only signiﬁcant sections of the graph should be taken in account and the regions beyond the cone of inﬂuence should not be explicated (Araghi *et al.* 2017). From the results, the WTC obtained for the NDVI–LST revealed some important in-phase patterns and the LST led the NDVI at 8–16 months band. Also, anti-phase patterns could be seen at 16–32 months band, and the NDVI led the LST in this band. At 8–16 months band, both in-phase and anti-phase relationships could be observed. The WTC of the NDVI–TCI showed anti-phase patterns at 8–16 and 16–32 months bands, and the LST led the NDVI. Some small areas with anti-phase relationships were detected at 4–8 months bands, in which the NDVI led the LST. The WTC of the SPI–VCI showed that the arrows were in-phase in most cases and there were positive relationships between the SPI and the VCI. The SPI–TCI wavelet analysis showed that in the 8–16 months band, the SPI led the TCI in all stations. For the 4–8 months band, in the S10 and S13 stations, there were in-phase relationships of the SPI and the TCI; while, for the S12 and S21 stations, an anti-phase relationship was found. Also, at the 32–64 months band, in-phase patterns could be seen for S10, S12, and S13, in which the SPI led the TCI. For the SPI–VHI, some areas with in-phase relationships were observed at the 4–8, 8–16, and 32–64 months bands. At the 8–16 months band, the VHI led the SPI and at the 32–64 months band, the SPI led the VHI. In general, it can be stated that the applied satellite-based drought indices have significant effects on the SPI drought index in the study area.

## CONCLUSIONS

Drought as one of the severe and frequently occurring natural hazards has negative impacts on human life. In this study, the spatiotemporal pattern of drought indices based on *in situ* observations and satellite datasets was investigated for the selected area. For meeting the objectives of this study, 29 *in situ* gauges as well as CMAP and MODIS satellite products were used. For having a correct vision of time-series decomposition via the MODWT, different types of wavelet filters were tested, and an optimal level for time-series decomposition was selected. The K-means and SOM clustering techniques were applied to regionalize the selected area. It was shown that the method based on the combined MODWT–EEMD–Energy–K-means technique is a robust approach for drought regionalization. The results showed that the stations located in a specific cluster were spread in different areas and the proximity of the stations to each other was not the sole basis of clustering. The satellite-derived indices including the NDVI, LST, VCI, TCI, and VHI time series were applied in order to investigate drought conditions and their inter-correlations. The results showed a negative relationship between the LST and the NDVI and between the TCI and the VCI and a positive relationship between the VHI and SPI time series. From the VHI time series point of view, the moderate and mild droughts occurred in the region during the selected time period. In the final step, the relationships between SPIs of the central stations of clusters and the NDVI and LST products were assessed via the WTC method. Based on the results obtained for the WTC analysis, the satellite-based indices have significant statistical correlations with the ground-based drought index within a certain period. The proposed method provides valuable references for drought assessment and management, specifically at the regional level. However, future work in the area can be done using a larger dataset that might yield better results. For the purpose of having more detailed and thorough investigations, it is also suggested to take other variables into consideration; for instance, soil moisture, other remotely sensed parameters, and large-scale processes. At some regions, the parameters that are correlated might be different based on the general precipitation cycle. Certain weather phenomenon such as the El Nino may also need to be taken into consideration.

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.