Predicting ecosystem shift in a Salt Lake by using remote sensing indicators and spatial statistics methods (case study: Lake Urmia basin)

Nadia Abbaszadeh Tehrani; Milad Janalipour

doi:10.4491/eer.2020.225

Environ Eng Res > Volume 26(4); 2021 > Article

Tehrani and Janalipour: Predicting ecosystem shift in a Salt Lake by using remote sensing indicators and spatial statistics methods (case study: Lake Urmia basin)

Research

Environmental Engineering Research 2021; 26(4): 200225.

Published online: August 9, 2021

DOI: https://doi.org/10.4491/eer.2020.225

Predicting ecosystem shift in a Salt Lake by using remote sensing indicators and spatial statistics methods (case study: Lake Urmia basin)

Nadia Abbaszadeh Tehrani^†

, Milad Janalipour

Department of Aerospace Application for Environment, Aerospace Research Institute, Ministry of Science, Research and Technology, Tehran, IRAN

^†Corresponding author: Email: Tehrani@ari.ac.ir, Tel: +982188366030, Fax: +982188366030

Received May 9, 2020 Accepted August 5, 2020

(open-access):

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

The consequences of unsustainable human activities on the environment are often delayed, when it is too late to compensate. New approaches are based on the use of “spatial statistics” of leading indicators to measure the “critical slowing down” in a degraded ecosystem, when it is reaching to a tipping point. This research predicts the tipping points in the ecosystem of Lake Urmia Basin (LUB) based on spatial statistics. By Remote Sensing (RS) indicators, their effectiveness in assessing the state of the ecosystem was evaluated in a 16-years period (2002–2017). Seven spectral indicators (NDVI, NDWIv, NDWIw, NDSI, SRDI, NMDI and MVWR) were extracted from ten MODIS images. Ability of the indicators to identify critical point in time-series was investigated by five spatial statistic methods (Moran’s-I, Getis-Ord-Gi, Geary’s-C, variance, and skewness). The results showed that Moran’s-I is more successful in predicting the ecosystem tipping point(s) in comparison with other methods. In addition, the ability to predict ecosystem trends by the autocorrelation of MVWR is higher than other indicators. According to results, the tipping points of LUB occurred in the years of 2008 to 2010 and 2015. For further studies, it is recommended to use radar indicators for identifying tipping points of the similar vulnerable ecosystems.

Keywords: Early warning signals, Ecosystem shift, MODIS, Remote sensing, Spatial indicator, Urmia Salt Lake

1. Introduction

In recent decades, most of the natural ecosystems of the world are at risk due to the human destructive activities. The impacts are often detectable when the ecosystem has passed from its acceptable thresholds and often entered into irreversible conditions or shifted to another stable state. According to [1] “Abrupt shift between ecological states may lead to ecological or economic losses. This happen when ecosystems reach to a tipping point, at which they may rapidly recognize into an alternative state with contrasting features”. These shifts could be predict in advance by leading indicators for critical transitions [2, 3]. It seems that remote sensing (RS) and Geographical Information System (GIS) can play an important role in predicting the ecosystem status before entering to the irreversible stage. “Alternative stable state theory” suggests that a change in ecosystem conditions can result in an abrupt shift in the state of the ecosystem. Ecosystems are often resistant to state shifts, this behavior is called “ ecosystem resilience” [4].

Monitoring the leading indicators of ecosystem shift in a vulnerable ecosystem can be effective in predicting ecosystem trends. This is one of the most important steps in natural resources and environmental management process. Several laboratory and field studies suggest that spatial statistics (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”.

When the resilience of an ecosystem decreases, it will recover slower from small disturbances until it is so close to a critical threshold that is also called a ‘tipping point’. In this situation, even a small stress can invoke a critical transition to another stable state where returning to the previous state can be difficult and sometimes impossible [6]. Several methods have been developed to measure this critical slowing down. These methods are based on spatial statistics of the leading indicators’ time-series.

Several studies have previously shown that spatial autocorrelation and standard deviation increase before a critical transition [3, 7–9]. When a stress drives the ecosystem towards the tipping point, the distribution of time-series become asymmetric, therefore, they show increased skewness. Ecosystems, which are close to the critical shift, have status variables that are timely and spatially auto-correlated. The increasing spatial coherence can be quantified by the spatial correlation function, between ecological states separated by a certain distance [1]. Increased recovery time may lead to stronger fluctuations around the equilibrium state and can increase spatial variance of the system [10].

Alibakhshini et al. [6] sign a critical change warning in a wetland ecosystem. In this research, the irreversible response of wetland ecosystem to pressure was monitored by RS and time-series data. They investigated the use of a series of RS data to track the key changes in the “Dorghe-Sangni” area near the Lake Urmia in Iran, which has undergone severe and destructive changes in recent decades. The used indicators include Normalized Difference Vegetation Index (NDVI), Modified Normalized Difference Water Index (MNDWI), Modified Vegetation/Water Ratio (MVWR) derived from the MODIS data series. They showed that autocorrelation in time lag1 of MVWR indicator can effectively indicate the early warning of critical changes in wetland ecosystem and the results are more acceptable compared to the use of NDVI and MNDWI. Seekell et al. (2012) [11] indicated that statistical anomalies (such as increased autocorrelation and variance) occur in time series spatial data when an ecosystem approaches to its tipping point. Statistical characteristics of autocorrelation and variance have been very effective in predicting ecosystem shift and can track changes that are very essential for adaptive management of ecosystems, biodiversity conservation and ecosystem services [12]. Previous studies on real ecosystems have shown that some indicators, especially the auto-correlation and variance of MVWR, can well predict critical changes in an heterogeneous ecosystem [6, 13, 14]. Due to unsustainable drainage of water resources supplying the Lake Urmia, serious environmental stresses have been induced to the lake that have threatened its health and sustainability over the past 15 years [15, 16]. Accurate and comprehensive data on the critical changes accrued in the lake can be obtained by RS. In this research, the validity of the selected remotely sensed spatial indicators of the ecosystem status was examined in Lake Urmia- as a real degraded aquatic ecosystem in a time period of 16 years (2002–2017). Accordingly, leading indicators were selected and developed by using appropriate spatial statistics and a method was introduced to monitor the ecosystem trend.

2. Material and Methods

2.1. Study Area

Lake Urmia, is a Salt Lake in northwest of Iran located between West and East Azerbaijan provinces (Fig. 1). Lake Urmia previously was the largest lake in the Middle East and the second largest Salt Lake on Earth [17, 18]. It is 4,183 feet above sea level, with geographic position between 45°E–46°E and 37°4′N–38°17′N, the lake has shrunk to 10% of its former size due to dam constructions and overuse of ground waters.

The current condition of Lake Urmia is the consequence of unsustainable development in the basin area for several decades. Fig. S1 shows the changes in the Lake Urmia’s surface area during the years 1976 to 2017 (Fig. S1 or Table S1 are presented in the supplementary file).

Water level of Lake Urmia has showed a sharp decrease in the years of 2002, 2008 and after 2010 over the past 30 years. In 2015, part of the lake has been restored with compensatory measures, but again, in 2017, the decreasing trend was observed. Various studies have shown that the “ecological surface level” of Lake Urmia is 1274.10 meters [19, 20]. Fig S2 shows water level changes of Lake Urmia from 1966 to 2010.

2.2. Data Used

Previous studies have used RS and interpretation of satellite imagery especially MODIS (Moderate Resolution Imaging Spectroradiometer) images to monitor ecosystem status and its resilience after facing a stress [21–23]. Here, Images of a particular month of the years (July) were imployed. Because the maximum changes of the indicators have been seen in the relevant month. Table S1 shows the selected MODIS images of Lake Urmia basin (LUB) from 2002 to 2017.

2.3. Methodology

In order to determine the effectiveness of selected indicators in predicting the changes in ecosystem status, their trends in time series were monitored by using spatial statistics. Analyzing the correlation between the selected indices and the degradation process in a real extra-saline aquatic ecosystem could help to choose the final spatial indicators of ecosystem shift. At first, the average of selected indices in LUB includes NDVI, SRDI, NDSI, NDWIW, NDWIV, NMDI and MVWR were calculated in time series to see if their average change can be considered as the early warning signals of ecosystem shift. Since the selected indicators had acceptable results, we selected them in this study [24].

The proposed method is summarized in four main steps:

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.

As the whole basin interactions affect Lake Urmia ecosystem and vice versa, it was chosen as the study area. In this study, indicators were selected based on the characteristics of the study area. Vegetation as well as water indicators were employed, since they are the main land cover classes in the study area. In addition, drought and salinity indices have been applied due to the dryness and salinity of the study area. For the first time, in this study, the probability of approaching to the “tipping point” in a Salt-Lake ecosystem is predicted by using remotely sensed indicators and their spatial statistics which are presented in Table 1.

2.3.1. Using spatial statistics

To predict ecosystem state, the recent approaches used spatial statistics. In previous studies, [6] have shown that metric-based methods (such as: Autocorrelation [29, 30],Variance [31] and Skewness [2]) applied on a time-series data of ecosystem status variable can precisely predict a critical catastrophic change in an ecosystem. For this reason, metric-based methods were employed to investigate the status of Urmia Lake ecosystem health.

2.3.1.1. Autocorrelation as a slowing down-based indicator

Increasing spatial integrity in an ecosystem can be quantified by the spatial correlation function of “Moran’s I” between the ecological states separated by a given distance (time). Autocorrelation can be calculated considering a certain time lag. Most of the studies on leading indicators have considered time lags of one time-step, or lag-1 autocorrelation [6]. In this study, Five methods including Variance, Skewness, Morans’I, Getis Ord Gi and Geary’s C were applied to estimate spatial autocorrelation.

a. Moran’s I method

Huge changes in an ecosystem may cause spatial variations. Moran’s I is a powerful method to measure spatial variations and is calculated according to Eq. (1) [32, 33].

(1)

I = \frac{n}{S_{0}} \frac{Σ_{i = 1}^{n} Σ_{i = 1}^{n} w_{i j} (x_{i} - \bar{x}) (x_{j} - \bar{x})}{Σ_{i = 1}^{n} {(x_{i} - \bar{x})}^{2}}

where x_i is the spectral value of pixel i obtained from spectral index and χ̄ is the average of the spectral value in a local area or entire of image. w_ij is the spatial weight between pixels i and j. Moreover, n is the total of pixels, I is Moran’s I and S₀ is obtained from Eq. (2).

(2)

S_{0} = Σ_{i = 1}^{n} Σ_{j = 1}^{n} w_{i j}

b. Getis-Ord Gi method

Getis-Ord Gi statistic is another method to compute spatial autocorrelation. Similar to Moran’s I index, Getis-Ord Gi (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.

(3)

G i = \frac{Σ_{j = 1}^{n} w_{i j} x_{j} - \bar{X} Σ_{j = 1}^{n} w_{i j}}{S \sqrt{\frac{n Σ_{j = 1}^{n} w_{i j}^{2} - (Σ_{j = 1}^{n} w_{i j})^2}{n - 1}}}

where X̄ and S are estimated from Eq. (4) and Eq. (5).

(4)

\bar{X} = \frac{Σ_{j = 1}^{n} x_{j}}{n}

(5)

S = \sqrt{\frac{Σ_{j = 1}^{n} x_{j}^{2}}{n} - {(\bar{X})}^{2}}

c. Geary’s C method

Geary’s C (C) method is used to compute spatial autocorrelation. It is sensitive the form of the graph of neighbor relationships According to equation, the Geary’s C method is computed according to Eq. (6) [32, 33].

(6)

C = \frac{(n - 1) Σ_{i} Σ_{j} w_{i j} {(x_{i} - x_{j})}^{2}}{2 S_{0} Σ_{i} {(x_{i} - \bar{x})}^{2}}

d. Variance and Skewness

Increasing the ecosystem recovery time leads to the more robust fluctuations around ecosystem balance. This will lead to the spatial variance before the critical transition. In such cases, we can observe a rise in skewness. Skewness indicates the extent to which the probability distribution is asymmetric. Reducing the recovery rate of a status variable can lead to a shift from equilibrium state and increases the variance. The skewness of the distribution of states is expected to increase not only if the system approaches a catastrophic bifurcation, but also if the system is driven closer to the basin boundary by an increasing amplitude of perturbation [2, 34]. The variance (SD) and skewness (SK) equations are calculated according to Eq. (7) and Eq. (8).

(7)

S D = \sqrt{\frac{1}{n - 1} Σ_{i = 1}^{n} {(x_{i} - \bar{x})}^{2}}

(8)

S K = \frac{\frac{1}{n} Σ_{i = 1}^{n} {(x_{i} - \bar{x})}^{3}}{\sqrt{\frac{1}{n} Σ_{i = 1}^{n} {(x_{i} - \bar{x})}^{2}}}

3. Results and Discussion

In order to test homogeneity of data, Mann-Kendall test was employed [35–37]. Z score obtained from the Mann-Kendall test showed changes of data. If Z score is higher than 1.96 or lower than −1.96, there is a great change in time-series. According to Table 2, a great change is observed in Moran index achieved from indicators. Moreover, it seems a great change is existing in MVMR, NDVI, and MDWIw, since main indicators have high Z score values. It must be considered that great changes cannot show tipping point and an indicator with high Z score value may be unable to find tipping point. For this reason, a comprehensive analysis about ability of indicators to find tipping point is presented in the following.

Fig. 3 shows the trend of changes of the selected indices including NDVI, SRDI, NDSI, NDWIw, NDWIv, NMDI and MVWR in LUB. Graphs indicate that there was a sudden change in 2008 in all of the selected indices. For all of the seven selected indices, the relevant spatial statistics were computed for 10 images in the period of 2002–2017. Due to the large volume of data and information, only the MVWR’s full results are presented here.

NDVI index was used to monitor vegetation cover changes. Fig. S3 shows the result of NDVI spatial statistics. The results showed that NDVI variance has a gradual decrease and increase respectively before and after the tipping point (2009 and 2015). It means that the minimum of NDVI was matched with the tipping point. A number of studies have shown that sometimes variance can be reduced when an ecosystem reaches to a tipping point [13]. Different studies have shown a continuous increase in autocorrelation and decreased/increased in variance [13]. In the case of Lake Urmia, the cause of variance reduction was the reduced ecosystem susceptibility and its failure to respond to the environmental fluctuations when approaching to the breakpoint. In the analysis of NDVI skewness curve, maximum increase occurred at tipping points. The reason for increased NDVI skewness in critical years in Lake Urmia may be the increased fluctuations around the state of equilibrium and irregularities in the changes.

NDSI is applied to monitor salinity, which ranges from −1 to +1. The higher NDSI value (the brighter regions in the image) indicates the more salinity. The areas around the boundary of Lake Urmia showed higher values over time, where were converted from water to saline area. Fig. S4 shows the result of NDSI spatial statistics.

The variance of NDSI has gradually declined over time as it reaches to the critical points (contrary to the expected increase in other studies). In other words, some ecosystems transform to several small areas after tension, and others are shrunk to one small area. The ecosystem of Lake Urmia seems to be closer to the second case, after perturbations.

Studying the time series of NDSI skewness showed a gradual increase and decrease respectively before and after the tipping point. The reason is that asymmetric changes and fluctuations around the mean of the indicator of salinity are seen due to a tension. In the review of NDSI autocorrelation in time-series the spatial autocorrelation (Moran’s I) has increased from at least five years before the tipping point. The reason could be the steady state of salinity condition which means the ecosystem’s failure to recover and reduce salt marshes between the two observations. It can be seen that among the three methods used to derive spatial autocorrelation, Moran’s I method gives better results in predicting the tipping point(s).

SRDI was developed to monitor the drought phenomenon. In SRDI, brighter gray grades for areas most affected by drought, and darker areas exhibit less variation in drought. Studying the extracted SRDI maps showed that the indicator amount was gradually increasing in the eastern and southern parts of the Lake. This is quite logical due to the conversion of water to bare land and salt marshes. The results of applying spatial statistics (variance, skewness, and autocorrelation) on the SRDI indicators are presented in Fig. S5.

SRDI variance gradually decreases over time as it approaches to the critical points, similar to the trends of the variance of other indicators. Studying the skewness of time-series of SRDI reveals that before the tipping point, the gradual decrease of the indicator (to the more negative values and in fact increase in its absolute value) was seen, and at the tipping points (especially at 2008; due to the global drought phenomenon at 2008 and its enormous regional impact on the ecosystem under study), the minimum of SRDI skewness has accrued. SRDI spatial autocorrelation (Moran’s I) was increasing during the several years before the breakpoint and among the three methods, was the best way to predict tipping point by using SRDI indicator.

NDWIw has been developed to monitor water bodies. By using this indicator, the water content of aquatic ecosystems such as lake and wetlands could be estimated. The range of the index is from −1 to +1 and the higher value indicates higher water content. The positive values related to the water zones and the zero or negative values related to vegetation or soil. Studying the time series of NDWIw variance showed that there is a gradual decrease and increase respectively, before and after the tipping point in the years 2008 and 2015. NDWIw skewness showed a big change in 2008, and the peak of that accrued in 2011. Fig. S6 shows the spatial statistics (autocorrelation, variance, and skewness) applied on NDWIw indicator.

NDWIv is used to monitor the water stress in the vegetation. The higher value of NDWIv shows the greater water content of the vegetation. In the extracted images, water and vegetation zones are seen more brightly. The results of applying spatial statistics on the NDWIv indicator are presented in Fig. S7.

Results showed that the NDWIv peak happened in 2010 and 2012. But, the variance of the indicator was decreasing as it approached to the tipping points-like the trend have seen in other indicators. Minimum of variance curve could be seen in the years 2008 and 2015 (breakpoints). NDWIv skewness chart indicates that the peak exactly matches the breakpoints (2010 and 2015). Among the three methods used to calculate autocorrelation, the best results were given by Moran’s I. It has been steadily raised as it approached to the tipping points, with peaks at critical years (2008 and 2013).

NMDI is a drought monitoring indicator. The results of this study on the relationship between drought and NMDI index are shown in Table S2. As can be seen, there is a relationship between drought and NMDI indicator in low vegetation conditions with leaf area index (LAI) less than 1. If the vegetation is dense (LAI = 2), this index will behave like NDWI and NDVI. The results of spatial statistics (variance, skewness, and autocorrelation) extracted from NMDI indicator are presented in Fig. S8.

The variance of NMDI was decreasing, reaching its lowest level in 2008, and then remaining almost constant after the critical shift of the ecosystem and the massive drought of 2008. Skewness, by 2011, had an upward trend toward negative numbers. The minimum has occurred in 2011 and was negative. It has increased since then, with the highest level in 2015. Autocorrelation towards the tipping points was increasing and the first peak point was seen in 2010. The second peak was observed in 2013 and third one in 2015 which is matched with the ecosystem’s tipping point.

MVWR was recommended to simultaneously consider vegetation and water condition. Due to the high capability of the indicator to show the status of heterogeneous ecosystems in earlier researches [6], it has also been applied in this study. In order to extract MVWR, at first, NDVI and MNDWI indicators have been extracted from the MODIS images. MVWR distribution in LUB during the years (2002–2017) has showed in (Fig. 4).

MVWR value varies from positive to negative. The higher and more positive values indicate higher levels of green vegetation and lower value, indicating more water content. Spatial statistics have been applied to MVWR and the results are presented in Fig. 5. The results show that the variance of the indicator has decreased from several years before the tipping point in 2008 (such as the pattern found in other indicators). But the value of its skewness has been steadily increased. After the tipping point, as the ecosystem state became stable, skewness was reduced. The results of the MVWR time-series skewness, can serve as a good benchmark for predicting the ecosystem’s critical shifts, due to showing their peak at the tipping point of 2008–2010 and 2015. The autocorrelation rate has an increasing trend until the breaking point of 2008–2010, 2012, and 2015, and can well predict the ecosystem behavior before reaching to the tipping points.

In order to test the effect of window size on tracking changes in the selected indicators, different window sizes have been used to calculate the variance of MVWR indicator. Fig. S9 compares the effect of the different window sizes on variations of the MVWR index in the time series.

In order to select the best window size, the standard deviation of MVWR variance was calculated in different years (Fig S10). The results showed that textural features extracted from a larger the size of the window are most appropriate to detect the tipping point. Although window size of 23 × 23 and 25 × 25 seem to be more suitable for tracking the changes by the variance of MVWR indicator, it was also observed that all windows follow a similar pattern.

Since 2000, the surface level of Lake Urmia has fallen below its ecological level and has continued its downward trend. Rate of the changes has increased since 2008 and the lake has reached to its minimum volume in 2015. Again with a restoration program in 2016, the lake condition was improved. The decline of the Lake surface area has started again since 2017. Since 2000 year, correlation of variance, skewness and autocorrelation of selected indicators with surface area changes in Lake Urmia are shown in Table S3 and Fig. S11.

The variance of different indicators had a high positive correlation with the surface area changes of the Lake. Meanwhile, the NMDI has the most correlation with Lake’s changes. The correlations of skewness for SRDI, NDSI, NMDI and MVWR indicators and NDVI, NDWIw, and NDWIv ones are positive and negative, respectively. The correlation between autocorrelation of indicators was generally negative and indicates the increase of autocorrelation by decreasing the surface area of the lake. Table S4 shows the correlation coefficient of variance, skewness and autocorrelation of selected indices with the surface level of Lake Urmia.

Fig. S12 compares the correlation of variance, skewness and autocorrelation with the Lake Urmia level. As can be seen, a high positive correlation is observed between the variance of most indices and the surface level of lake. Autocorrelation of indicators often had a moderate negative correlation with the level of lake. NMDI and NDSI had the highest value, while NDWIv and NDWIw had the lowest one. The correlation of skewness varied with the level changes from negative to positive, the highest and smallest values belong to NDVI and NDWIw, and NDWIv, respectively. The results indicate the appropriate correlation between variance and autocorrelation of the selected indicators with the surface level and these leading indicators could be used confidently in the Lake changes analysis.

In order to prove the effectiveness of the selected indicators in predicting sudden changes in the ecosystem situation, it is necessary to assess time correlation of each indicator in a period. For this purpose, the correlation coefficient of each indicator with one-time interval (time delay 1) was extracted. Fig S13 shows the correlation coefficient of each indicator with a subsequent time interval (time delay1 (. Results show that the minimum values of the NDVI, SRDI, NDWIv, and MVWR time correlation curves are occurred in 2010–2008 and 2015–2013. Therefore, it can be deduced the mentioned indicators can predict the tipping points. Moreover, the minimum values of the NDSI, NMDI, and NDWIw time correlation curves haven not matched with the tipping points.

In Table S5, the behavior of variance of the selected indicators regarding the tipping point is compared. Trend of Variance towards the tipping point is decreasing. The minimum point of variance curve in the water indicators such as NDWIw and NDWIv could be used to predict the tipping point. Peak of variance of other indices such as NDVI, NDSI and SRDI happened in 2010 a year after the tipping point and therefore cannot be used to predict the tipping point.

In Table S6, trend of skewness of the selected indicators towards the tipping point(s) is analyzed. The skewness of NDVI can identify the tipping point. Skewness of NDWIv and NDSI has been rising steadily up to 2010, and then decreased and could be successful in predicting tipping points. The highest value of skewness of NMDI is in 2011, and then decreases. Hence, it is successful in predicting the trend, but its maximum value is in 2011 instead of 2008. Therefore, it cannot correctly predict the critical point. MVWR’s skewness will also increase by 2010, reaching its lowest level in 2013. In general, skewness was successful in predicting breakpoints in some of the indicators (NDWIv and NDSI), while it wasn’t in some others (NMDI, MVWR, NMDI and NDVI).

In Table 3, the maximum and minimum values of the autocorrelation curves of selected indices of LUB have been compared. The autocorrelation of NDVI is able to identify the initial tipping point (in 2008). SRDI’s autocorrelation has been steadily rising from several years before the tipping point (from 2002 to 2010), and then decreased. Then, it could be successful in predicting the ecosystem trend. NDSI’s autocorrelation has increased from 2002 to 2010, but it could predict the general trend of the ecosystem. Investigating the autocorrelation curves of NDWIw and NDWIv shows that the peak of these indicators happened in the years of 2008 and 2013 that can be used to determine the ecosystem status. NMDI and MVWR have shown an increasing trend since 2006, reaching their peak in 2010, 2013 and 2015. Due to the continuous increase before the breakdown, they can be considered as leading indicators for predicting changes in ecosystem status.

4. Conclusions

In this study, spatial statistics methods (variance, skewness, and autocorrelation) were applies on seven selected remote sensing indicators, including: vegetation (NDVI), drought (NMDI), salinity (NDSI), water (NDWIv, NDWIw) and Vegetation/water (MVWR) indices to study the trend of their changes in a 16-years period (2002–2017) in LUB. It was proved that applying the spatial autocorrelation method (Moran’s I) has better results in comparison with two other methods (spatial variance and spatial skewness) in predicting the ecosystem tipping point. In addition, autocorrelation of MVWR indicator provides a better prediction of the tipping point/s in the LUB compared to other indicators.

The results indicated that in the mixed and heterogeneous ecosystem of the LUB, with different land cover, water and vegetation zones, among the selected indicators, the MVWR index, which is a mixed index of water (NMDWw) and plant (NDVI) indicator, provides a better result in predicting the ecosystem tipping point/s. It was revealed that the spatial statistics (Moran’s I) of seven selected indicators have increased, until reaching to the tipping point. Increasing the autocorrelation of an indicator in a time series, can indicate that the ecosystem is close to its tipping point. The achievements were consistent with the results of previous studies as well. According to each ecosystem conditions, it is necessary to introduce the indicators of change prediction. In this study, due to the unique conditions of LUB, other indicators such as NDSI and SRDI have also been used to evaluate their ability to predict the state of the ecosystem. The results showed the good performance of these indicators, especially their autocorrelation in prediction the critical points. In prioritizing the indicators to predict the tipping point/s, spatial autocorrelation, spatial skewness and spatial variance have better performance, respectively. Based on the accuracy of the results, the following indicators are recommended to predict ecosystem shift in the study area, respectively: MVWR, NMDI, NDWIv, NDWIw, NDSI, SRDI, NDVI. Therefore, the vegetation /water compound indicator, as well as drought and even salinity indices, were better than NDVI, in predicting tipping points in the LUB. In general, it seems that it is possible to monitor ecosystem health using RS based indicators and in order to predict the condition of any ecosystem, according to its nature and the type of stresses imposed on it, it is necessary to select appropriate and effective indicators that consider all components involved in the health of the ecosystem in a comprehensive and integrated manner. LUB is an ecosystem containing a unique Salt Lake that the impacts of climate change and drought phenomenon combined with the stresses resulting from unsustainable human activities in the form of land use development, agriculture, dam construction, water drainage, endangering its health and accelerating its regression process. In this study, in addition to common land surveying indicators such as NDVI, indicators that show the health of the entire catchment area have been introduced. For further study, climatic and complementary data, such as land use changes can be used to monitor the trend of changes in indicators. The present research was conducted on a Salt Lake, which had already experienced the tipping points. In future studies, the proposed method can be employed in other ecosystems. Finally, it is recommended to use radar remote sensing data and determine its appropriate indicators for identifying tipping points of the similar vulnerable ecosystems.

Supplementary Information

eer-2020-225-suppl.pdf

Acknowledgment

The authors would like to thank Iranian Space Agency (ISA) for providing MODIS data sets and Aerospace Research Institute (ARI) for its financial support for this study.

Notes

Author Contributions

Dr. N.A.T. (assistant professor) designed experiments and investigated outputs and wrote the manuscript. Dr. M.J. (assistant professor) designed experiments and applied the designed experiments and revised the manuscript.

References

1. Kefi S, Guttal V, Brock WA, et al. Early warning signals of ecological transitions: methods for spatial patterns. PloS one. 2014;9:e92097.

2. Scheffer M, Bascompte J, Brock WA, et al. Early-warning signals for critical transitions. Nature. 2009;461:53.

3. Dakos V, Carpenter SR, Brock WA, et al. Methods for detecting early warnings of critical transitions in time series illustrated using simulated ecological data. PloS one. 2012;7:e41010.

4. Holling CS. Resilience and stability of ecological systems. Annu Rev Ecol Syst. 1973;4:1–23.

5. Dakos V, van Nes EH, Donangelo R, Fort H, Scheffer M. Spatial correlation as leading indicator of catastrophic shifts. Theor Ecol. 2010;3:163–174.

6. Alibakhshi S, Groen TA, Rautiainen M, Naimi B. Remotely-sensed early warning signals of a critical transition in a wetland ecosystem. Remote Sens. 2017;9:352.

7. Lamsal K. Identifying potential critical transitions in a forest ecosystem using satellite data [Thesis]. Sweden: Lund University GEM thesis series; 2017.

8. Weaver I, Dyke J. Early warning signals in complex ecosystems. Earth Syst Dynam Discuss. 2015;6:2507–2542.

9. Lenton T, Livina V, Dakos V, Van Nes E, Scheffer M. Early warning of climate tipping points from critical slowing down: comparing methods to improve robustness. Phil Trans R Soc A. 2012;370:1185–1204.

10. Guttal V, Jayaprakash C. Spatial variance and spatial skewness: leading indicators of regime shifts in spatial ecological systems. Theor Ecol. 2009;2:3–12.

11. Seekell DA, Carpenter SR, Cline TJ, Pace ML. Conditional heteroskedasticity forecasts regime shift in a whole-ecosystem experiment. Ecosystems. 2012;15:741–747.

12. Burthe SJ, Henrys PA, Mackay EB, et al. Do early warning indicators consistently predict nonlinear change in long-term ecological data? J Appl Ecol. 2016;53:666–676.

13. Dakos V, Van Nes EH, D’Odorico P, Scheffer M. Robustness of variance and autocorrelation as indicators of critical slowing down. Ecology. 2012;93:264–271.

14. Carpenter S, Brock W, Cole J, Kitchell J, Pace M. Leading indicators of trophic cascades. Ecol Lett. 2008;11:128–138.

15. Hoseinpour M, Fakheri Fard A, Naghili R. Death of Urmia Lake, a silent disaster investigating causes, results and solutions of Urmia Lake drying. In : 1st International Applied Geological Congress, Department of Geology, Islamic Azad University-Mashad Branch; 26–28 April 2010; Iran. p. 700–704.

16. Ghale YAG, Altunkaynak A, Unal A. Investigation anthropogenic impacts and climate factors on drying up of Urmia Lake using water budget and drought analysis. Water Resour Manag. 2018;32:325–337.

17. Zarghami M. Effective watershed management; case study of Urmia Lake, Iran. Lake Reserv Manag. 2011;27:87–94.

18. Hassanzadeh E, Zarghami M, Hassanzadeh Y. Determining the main factors in declining the Urmia Lake level by using system dynamics modeling. Water Resour Manage. 2012;26:129–145.

19. Abbaspour M, Nazaridoust A. Determination of environmental water requirements of Lake Urmia, Iran: an ecological approach. Int J Environ Stud. 2007;64:161–169.

20. Dariane A, Eamen L. Finding the Causes and Evaluating Their Impacts on Urmia Lake Crisis Using a Comprehensive Water Resources Simulation Model. J Hydraul Struct. 2017;3:62–77.

21. Liu X, Jiang W, Li J, Wang W. Evaluation of the Vegetation Coverage Resilience in Areas Damaged by the Wenchuan Earthquake Based on MODIS-EVI Data. Sensors. 2017;17:259.

22. Zhang K. Regime shifts and resilience in China’s coastal ecosystems. Ambio. 2016;45:89–98.

23. Simoniello T, Lanfredi M, Liberti M, Coppola R, Macchiato M. Estimation of vegetation cover resilience from satellite time series. Hydrol Earth Syst Sci Discuss. 2008;5:511–546.

24. Babaei H, Janalipour M, Tehrani NA. A simple, robust, and automatic approach to extract water body from Landsat images (case study: Lake Urmia, Iran). J Water Clim Chang. 2019;

25. Gao BC. Normalized difference water index for remote sensing of vegetation liquid water from space. Imaging Spectrometry. 1995;2480:225–236.

26. Richardson AD, Duigan SP, Berlyn GP. An evaluation of noninvasive methods to estimate foliar chlorophyll content. New Phytol. 2002;153:185–194.

27. Iii JE, McLeod KW. Indications of relative drought stress in longleaf pine from thematic mapper data. Photogramm Eng Remote Sensing. 1999;65:495–501.

28. Wang L, Qu JJ. NMDI: A normalized multi-band drought index for monitoring soil and vegetation moisture with satellite remote sensing. Geophys Res Lett. 2007;34:

29. Wissel C. A universal law of the characteristic return time near thresholds. Oecologia. 1984;65:101–107.

30. Strogatz SH. Nonlinear Dynamics and Chaos with Applications to Physics, Biology, Chemistry, and Engineering. 1st edNew York: Perseus BooksPublishing; 1994.

31. Carpenter SR, Brock WA. Rising variance: a leading indicator of ecological transition. Ecol Lett. 2006;9:311–318.

32. Getis A, Ord JK. The analysis of spatial association by use of distance statistics. Perspectives on Spatial Data Analysis, ed. Berlin: Springer; 2010. p. 127–145.

33. Goodchild MF. Spatial autocorrelation. Geo Books, University of Western Ontario; Canada, Ontario: Canada: 1986. p. 56.

34. Dakos V, Scheffer M, van Nes EH, Brovkin V, Petoukhov V, Held H. Slowing down as an early warning signal for abrupt climate change. Proc Natl Acad Sci. 2008;105:14308–14312.

35. Kousari MR, Zarch MAA, Ahani H, Hakimelahi H. A survey of temporal and spatial reference crop evapotranspiration trends in Iran from 1960 to 2005. Clim Change. 2013;120:277–298.

36. Ahani H, Kherad M, Kousari MR, et al. An investigation of trends in precipitation volume for the last three decades in different regions of Fars province, Iran. Theor Appl Climatol. 2012;109:361–382.

37. Jaagus J. Climatic changes in Estonia during the second half of the 20th century in relationship with changes in large-scale atmospheric circulation. Theor Appl Climatol. 2006;83:77–83.

Fig. 1

Location of the study area and false RGB color image of MODIS in 2004.

Fig. 2

Research Flowchart based on RS indicators to find tipping point of Lake Urmia ecosystem.

Fig. 3

Selected indicators’ trend in LUB including (a) NDVI (b) SRDI (c) NDSI (d) NDWIw (e), NDWIv, (f) NMDI and (g) MVWR in time series (from year 2002 to 2017) (Horizontal and vertical axes represent years and value, respectively).

Fig. 4

MVWR map of LUB extracted from the time series of the MODIS images during the years 2002 to 2017.

Fig. 5

Spatial statistics (Autocorrelation, Variance, Skewness) applied to MVWR indicator. (Horizontal and vertical axes represent years and value, respectively).

Table 1

Selected Indicators to Predict Tilling Point in LUB

No	Abbreviation	Full Name	Formula	MODIS bands	Reference	Explanation
1	NDVI	Normalized Difference Vegetation Index	$N D V I = \frac{N I R - R e d}{N I R + R e d}$	$N D V I = \frac{b 2 - b 1}{b 2 + b 1}$	[6]	To find vegetation change of the area
2	NDWI-V	Modified Normalized Difference Water Index(leaves)	$N D W I v = \frac{N I R - S W I R}{N I R + S W I R}$	$N D W I v = \frac{b 2 - b 6}{b 2 + b 6}$	[25]	It is robust index to find water of vegetation
3	NDWI-W	Modified Normalized Difference Water Index	$N D W I w = \frac{G r e e n - N I R}{G r e e n + N I R}$	$N D W I w = \frac{B 4 - B 2}{B 4 + B 2}$	[25]	since a great change over water body is occurred, this index can be useful
4	NDSI	Modified Normalized Difference Salinity Index	$N D S I = \frac{S W I R 1 - S W I R 2}{S W I R 1 + S W I R 2}$	$N D S I = \frac{B_{6} - B_{7}}{B_{6} + B_{7}}$	[26]	Salinity index is helpful to find changes water to saline
5	SRDI	Simple Ratio Drought Index	$S R D I = \frac{M I R}{N I R}$	$S R D I = \frac{B 7}{B 5}$	[27]	It is a simple index, which uses two effective bands
6	NMDI	Normalized Multi-band Drought Index	$N M D I = \frac{N I R - (S W I R 1 - M I R)}{N I R + (S W I R 1 - M I R)}$	$NMDI = \frac{B 2 - (B 6 - B 7)}{B 2 + (b 6 - b 7)}$	[28]	It is a powerful index for drought monitoring
7	MVWR	Modified Vegetation/ Water Ratio	$M V W R = ln (\frac{N D V I + 1}{M N D W I + 1})$	-	[6]	This index can consider water and vegetation simultaneously

where, NIR = Near Infrared band. Green= green band, Red = red band, MIR = Middle Infrared band. SWIR= Shortwave infrared band, SWIR1 = Band in range (1628–1700) nm, SWIR2 = Band in range (2105–2155) nm

Table 2

Z Score Obtained from Mann-Kendall Test to Show Homogeneity of Data

Indicator	Mean	Variance	Moran	Getis	Geary
MVMR	1.64	−2.97	2.26	−3.09	−0.04
NDSI	0.78	−0.28	−0.61	2.7	−2.26
NDVI	2.1	−3.52	3.50	−2.34	−2.10
NDWIv	2.75	−1.99	1.84	1.02	0.13
NDWIw	−2.96	−3.19	3.17	1.77	0.70
SRDI	−1.68	−3.17	2.59	3.66	1.48
NMDI	−2.22	1.80	−1.11	0.58	3.64

Table 3

Comparison of Maximum and Minimum of the Autocorrelation Curves of Selected Indices in LUB

Selected indicators	Trend of autocorrelation before and after the tipping point(s)	Minimum (year)	Maximum (year)	Tipping point of the ecosystem (year)
NDVI	Ascending-Descending	2004,2006	2008,2012	2008–2015
SRDI	Ascending-Descending	2002,2006,2013	2010,2012	2008–2015
NDSI	Ascending-Descending	2002,2017	2010	2008–2015
NDWIw	Ascending-Descending	2006	2008,2013	2008–2015
NDWIv	Ascending-Descending	2006	2008,2013	2008–2015
NMDI	Ascending-Descending	2004,2006	2010,2013	2008–2015
MVWR	Ascending-Descending	2004,2006	2010,2015	2008–2015