2.3. Methodology
EVI and LST images were used to calculate Vegetation Condition Index (VCI) and Temperature Condition Index (TCI), respectively, in order to prepare VHI. Then, VHI was employed to determine PVVI. It should be mentioned that EVI in each pixel has different values in different months. Thus, the maximum EVI, i.e. when the vegetation cover reaches its annual peak, should be identified and considered in each pixel. In the study area, the highest amount of EVI is in April and May, so the average EVI of these two months was considered as the maximum annual EVI. Moreover, the data of precipitation and temperature were used to calculate SPEI. The SPEI calculations were carried out and coded in MATLAB R2021b for each station on 3, 6, 9 and 12-month time scale. Then, they were interpolated across the study basin using Inverse Distance Weighted (IDW) model in ArcMap 10.8.1 to obtain SPEI maps. The pixel size and position were considered as EVI maps.
Then, Pearson correlation coefficient and linear regression was applied to investigate the relationship between EVI and multitemporal SPEI. The maximum value of the correlation for each pixel was considered to evaluate the response of vegetation cover to different SPEI time scales. The slope of the linear regression indicates the trend of the EVI and multitemporal SPEI (Shi et al. 2021).
Afterward, the predicted EVI was produced for each pixel based on the correlation between EVI, obtained from satellite images (observed EVI) and multitemporal SPEI. Then, the residual was obtained from the difference of predicted EVI and observed EVI. The effects of climate change and human activities on vegetation cover is distinguished (Li et al. 2012) and classified (Sun et al. 2015) using residual analysis according to Table 2.
Table 2
Methods for classification of climate change and human activities roles (Sun et al. 2015)
Vegetation trend
|
Slope(predicted)
|
Slope(residual)
|
Vegetation cover is induced by:
|
Increasing (slopeobserved > 0)
|
> 0
|
> 0
|
Both climate change and human activities
|
> 0
|
< 0
|
Climate change
|
< 0
|
> 0
|
Human activities
|
Decreasing slope(observed) < 0)
|
< 0
|
< 0
|
Both climate change and human activities
|
< 0
|
> 0
|
Climate change
|
> 0
|
< 0
|
Human activities
|
Finally, the probability of vegetation vulnerability was compared in where affected by the climate change or human activities. The Pearson correlation coefficient, linear regression and residual analysis were carried out in TerrSet 2020 v.19.0.6. Figure 3 illustrates the flowchart of the methodology in this study.
PVVI
EVI was used to obtain the VCI according to Eq. 1. VCI compares current vegetation cover with previous values (Rimkus et al. 2017) and helps assess environmental moisture deficit (Jain et al. 2010). Therefore, this index is one of the suitable methods for evaluating and analyzing drought using satellite images.
VCI =\(\frac{{\text{E}\text{V}\text{I}}_{\text{i}}-{\text{E}\text{V}\text{I}}_{\text{m}\text{i}\text{n}}}{{\text{E}\text{V}\text{I}}_{\text{m}\text{a}\text{x}}-{\text{E}\text{V}\text{I}}_{\text{m}\text{i}\text{n}}}\) × 100 (1)
Where; EVIi is the EVI value of a pixel in year i, and EVImax and EVImin are the maximum and minimum values of EVI value of the same pixel from 2001 to 2019, respectively.
In the next step, TCI is calculated according to Eq. 2. This index makes it possible to investigate the response of vegetation to thermal stress. This index is suggested by (Kogan 1995).
TCI = \(\frac{{\text{L}\text{S}\text{T}}_{i}-{\text{L}\text{S}\text{T}}_{\text{m}\text{i}\text{n}}}{{\text{L}\text{S}\text{T}}_{\text{m}\text{a}\text{x}} + {\text{L}\text{S}\text{T}}_{\text{m}\text{i}\text{n}}}\) × 100 (2)
Where; LST is Land Surface Temperature, LSTi is the LST value of a pixel in year i, and LSTmax and LSTmin are the maximum and minimum values of LST value of the same pixel from 2001 to 2019, respectively.
In next step, VHI is obtained using VCI and TCI. This index reflects the health of vegetation cover (Kundu et al. 2016) and is employed to diagnose the occurrence, severity and duration of drought (Karnieli et al. 2010). VHI considers the characteristics of the ecosystem in terms of fluctuations of maximum and minimum vegetation and land surface temperature (Bento et al. 2018). VHI is computed as Eq. 3.
VHI = \({\alpha }\) × VCI + (1 - \({\alpha }\)) × TCI (3)
Where; \(\alpha\) is the weight of TCI and VCI and its amount depends on the temperature and humidity conditions of the area (Tran et al. 2017), but as the needs of plants are not clear during the growth cycle, it is considered 0.5 (Cong et al. 2017). Finally, VHI values are classified according to Table 3.
Table 3
Classification of VHI (Kogan 1995; Le Hung & Hoai 2015) and the weight of classes for PVVI (Heydari Alamdarloo et al. 2018)
VHI Class
|
VHI Value
|
Class of drought intensities
|
Weight classes in PVVI
|
1
|
0 ≤ VHI < 10
|
Extreme Drought
|
4
|
2
|
10 ≤ VHI < 20
|
Severe Drought
|
3
|
3
|
20 ≤ VHI < 30
|
Moderate Drought
|
2
|
4
|
30 ≤ VHI < 40
|
Mid Drought
|
1
|
5
|
40 ≤ VHI
|
No Drought
|
0
|
After estimating and classifying VHI, in order to determine the probability of vegetation vulnerability, PVVI is estimated. Given that the occurrence of different drought intensities in the probability of vegetation vulnerability in an area are not the same value, for quantification and calculation of PVVI, each class of VHI is weighted (Table 3) and then PVVI is determined using Eq. 4:
PVVI= \(\sum _{i=1}^{5}{W}_{i}\) × \({P}_{Ci}\) (4)
Where; PVVI is Probability of Vegetation Vulnerability Index, Wi is weight of class i and PCi is probability of occurrence of class i which is calculated using Eq. 5:
\({\text{P}}_{\text{C}\text{i}}=\frac{{\text{n}}_{\text{C}\text{i}}}{\text{N}}\) × 100 (5)
Where; \({n}_{Ci}\) is the number of occurrences of class i in a pixel during 2001–2019 and N is total number of VHI (N = 19).
Depending on the weights of the classes, PVVI ranges from 0 and 400, and the higher the value of this index, the more likely the vegetation is vulnerable. In this paper, the map of this index was classified into 5 classes based on natural break method in ArcMap 10.8.1.
SPEI
The monthly SPEI was estimated to study the trend of climate changes using precipitation and potential evapotranspiration. This index shows the difference between precipitation and reference evapotranspiration, and thus reflects the water deficit, which is a vital factor for vegetation activity in this region (Shi et al. 2021; Vicente-Serrano et al. 2010). Therefore, in this research, SPEI was used to investigate the response of vegetation to climate changes and human activities.
To estimate potential evapotranspiration, Thornthwaite (1948) equation was used:
$$PET=16{N}_{m}{\left(\frac{10{T}_{m}}{I}\right)}^{a}$$
6
\({N}_{m}\) is the correction coefficient that is determined based on the month and the latitude of the studied area, \({T}_{m}\) is average temperature of the month in degrees Celsius, I is thermal index calculated for the whole year and α is a coefficient computed based on I.
Next, the monthly water deficit or surplus (Eq. 8) and after that, \({X}_{i,j}^{k}\) was calculated.
$${D}_{i}={P}_{i}-{PET}_{i}$$
7
$${X}_{i,j}^{k}=\left\{\begin{array}{c}\sum _{l=13-k+j}^{12}{D}_{i-1,l}+\sum _{l=1}^{j}{D}_{i,j} if j<k\\ \sum _{l=j-k+1}^{j}{D}_{i,j} if j\ge k\end{array}\right.$$
8
Where; i is the respective year, j is the respective month and k is time scale (k = 3, 6, 9 and 12).
In order to standardize the series of water deficit or surplus observations to determine SPEI according to Vicente-Serrano et al. (2013), three-factor log-logistic distribution was used as the most appropriate adaptive distribution based on Eq. (9):
$$f\left(x\right)=\frac{\beta }{\alpha }{\left(\frac{x-\gamma }{\alpha }\right)}^{\beta -1}{\left[1+\left(\frac{x-\gamma }{\alpha }\right)\right]}^{-2}$$
9
Where; α, β, and γ are scale, shape, and boundary factors, respectively, for D values in the interval γ < D < ∞. Also, Eq. (10) shows the probability distribution function of the observation series of water deficit or surplus according to the log-logistic distribution.
$$F\left(x\right)={\left[1+{\left(\frac{x-\gamma }{\alpha }\right)}^{\beta }\right]}^{-1}$$
10
Eventually, based on F(x), the SPEI was computed as Eq. (11):
$$SPEI=W-\frac{{C}_{0}+{C}_{1}W+{C}_{2}{W}^{2}}{1+{d}_{1}w+{d}_{2}{w}^{2}+{d}_{3}{w}^{3}}$$
11
Where; \({C}_{0}=2.515517\), \({C}_{1}=0.802853\), \({C}_{2}=0.010328\), \({d}_{1}=1.432788\), \({d}_{2}=0.18926\) and \({d}_{3}=0.001308\). \(W=\sqrt{-2{ln}P}\) if P=1-F(x)\(\le 0.5\) but if P > 0.5, P was replaced by 1-P, where; P represents precipitation.
SPEI has positive and negative values, the more negative this value indicates drought and the positive value of this index indicates wet condition.