### 1. Introduction

*i.e.*, variance, skewness, and autocorrelation) of leading indicators could predict tipping point prior to transitions [1]. Due to increased recovery time to local equilibrium after a stress, neighboring units become more like each other when a system approaches a bifurcation point,

*i.e.*, they become increasingly correlated [5].The theories of critical change suggest that leading indicators change predictably before a tipping point and thus could be determined as the early signals of critical change. These changes are called “critical slowing down”.

### 2. Material and Methods

### 2.1. Study Area

### 2.2. Data Used

### 2.3. Methodology

Appropriate MODIS images were selected. Bands 1 and 2 with spatial resolution of 250 meters were resampled to 500 meters in order to create an image with a size equal to size of five bands with spatial resolution of 500 meters and then all seven bands were stacked in one layer. Finally, the study area was masked using the boundary layer of the basin.

Seven spectral indicators were extracted from ten images ( features were produced). Then the textural and spatial correlation features were extracted from them, since they are powerful in this regard. Means of texture, spectral and spatial correlation features for each year was calculated [24].

The ability of the features to identify the ecosystem’s critical point was investigated by using three spatial statistic methods for each of seven leading indicators in time series (features were produced).

By reviewing valid ground data and existing records, validation of the results was carried out. Moreover, at this stage, the accuracy of the results was examined. The conceptual framework of the research is presented in Fig. 2.

#### 2.3.1. Using spatial statistics

### 2.3.1.1. Autocorrelation as a slowing down-based indicator

### a. Moran’s I method

##### (1)

$$I=\frac{n}{{S}_{0}}\frac{{\mathrm{\Sigma}}_{i=1}^{n}{\mathrm{\Sigma}}_{i=1}^{n}{w}_{ij}({x}_{i}-\overline{x})({x}_{j}-\overline{x})}{{\mathrm{\Sigma}}_{i=1}^{n}{({x}_{i}-\overline{x})}^{2}}$$*x*

*is the spectral value of pixel i obtained from spectral index and*

_{i}*χ̄*is the average of the spectral value in a local area or entire of image.

*w*

*is the spatial weight between pixels*

_{ij}*i*and

*j*. Moreover,

*n*is the total of pixels,

*I*is Moran’s I and

*S*

_{0}is obtained from Eq. (2).

### b. Getis-Ord Gi method

*Gi*) was computed based on the difference between pixel value and the mean of digital numbers (Eq. (3)) [32, 33]. It was frequently used in previous studies.

### c. Geary’s C method

### d. Variance and Skewness

*SD*) and skewness (

*SK*) equations are calculated according to Eq. (7) and Eq. (8).