**1. Introduction**

The spatial interface between two or more ecological regions and their material, energy, and structural and functional systems is called an ecotone or an ecological ecotone [1]. The desert-oasis ecotone is distributed along the periphery of an oasis and is characterized by a zone of desert vegetation that separates the extensive desert from the oasis [2]. The ecotone records the interaction and mutual transformation between the desert and oasis ecosystems [3] and serves as an ecological link connecting the two. A desert-oasis ecotone is a unique ecosystem between a desert and an oasis, usually characterized by low diversity, sparse cover, and dominance by perennial herbaceous grasses and semi-shrubs, such as *Phragmites australis*, *Tamarix ramosissima*, *Karelinia caspia*, and *Alhagi sparsifolia.* The ecotone can be used for ranching (of both livestock and wild animals); its vegetation can also increase the roughness of the underlying ground surface, thereby hindering the development of desertification and protecting the oasis from wind erosion and sand deposition [4–7]. At the same time, the desert-oasis ecotone is the interface between the oasis ecosystem and the desert ecosystem where energy, material, and information exchange occurs [8], which is highly sensitive to external environmental and human disturbances, affected easily by human activities, including the expansion of cultivated land

**Citation:** Sun, F.; Wang, Y.; Chen, Y.; Li, Y.; Zhang, Q.; Qin, J.; Kayumba, P.M. Historic and Simulated Desert-Oasis Ecotone Changes in the Arid Tarim River Basin, China. *Remote Sens.* **2021**, *13*, 647. https:// doi.org/10.3390/rs13040647

Academic Editor: Giovanni Chirico Received: 21 January 2021 Accepted: 6 February 2021 Published: 11 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Beijing 102206, China

and urbanization and precipitation [9]. A desert-oasis ecotone being a natural ecological barrier that prevents the desert from expanding into the oasis, its analysis provides an important indicator and an early warning of ecological changes. The desert-oasis ecotone also plays a large role in the development of the oasis economy. Therefore, the ecotone has been a topic of significant research in recent years. Previous research on this ecotone primarily focuses on characterizing and classifying the vegetation diversity, microclimate, and soil moisture, among other parameters in the ecotone and/or oasis [8,10–12]. The comparison of these studies has partially revealed the causes of formation and the ecological environment of ecotones in deserts and oases of different arid areas, thus providing the empirical basis for the ecological protection of desert-oasis ecotones.

The observed changes and development of ecotones are closely related to the dynamic evolution of the overall environment and climate of the region, as well as human influencing factors. Maintaining the stability and development of oases requires a comprehensive analysis of long-term trends and of the causes of the changes in the ecotone as well as within the entire basin. However, such studies are generally lacking. Previous studies have also shown that changes in the ecotone are closely related to changes in land-use type, which provides research ideas for predicting future changes in the ecotone. Thus, the future dynamic evolution of the ecotone can be inferred by predicting land-use changes.

There are many methods of simulating and predicting the evolution of land-use patterns, such as system dynamics, CLUE-S, artificial neural network (ANN), and cellular automata–Markov chain (CA-Markov) methods [13–16]. The system dynamics model is based on cybernetics, information theory, and system theory to analyze the drivers of land-use change. At present, the system dynamics simulation software STELLA has not been fully combined with the spatial analysis function of GIS to implement land-use change simulation and has not played its role as a powerful dynamics system. [17]. The CLUE-S model has limitations, and the setting of some parameters in the model mainly relies on expert knowledge, which will bring a high degree of subjectivity [13]. The ANN simulation of land-use change requires long simulation times, and the method does not provide the user with a specific evolution formula and contains large errors [18]. The CA-Markov model not only retains the advantages of the Markov model for long-term prediction but also integrates the ability of the cellular automata (CA) model to simulate complex spatiotemporal system changes. Thus, the CA-Markov model can better simulate land-use changes in time and space and has been widely used [19–21].

In the Tarim River Basin, artificial oasis [22] and desertification processes are increasing [23]. As a result, the desert-oasis ecotone has rapidly decreased in size and ecological concerns have become increasingly prominent. Meanwhile, the rapid advancement of urbanization in the Tarim River Basin, as well as the continuous development of the social economy, has led to significant changes in the pattern of land use, which have produced a series of impacts on the ecological environment of the ecotone.

The main purposes of this study are (1) to analyze the spatial and temporal variability and driving forces of the desert-oasis ecotone in the Tarim River Basin from 1990 to 2015, (2) to evaluate the applicability of the CA-Markov model, and (3) to further predict the near-future land-use changes in the Tarim River Basin. This study will deepen our knowledge of the evolution of the desert-oasis ecotone, which has important implications for the protection of the ecological environment in the arid zone and the construction of an ecological civilization in the Silk Road Economic Belt.

#### **2. Study Area**

The Tarim River Basin is an inland basin located far from the ocean in northwest China (Figure 1). The area is characterized by a temperate arid continental climate with scarce precipitation and strong evaporation. In the study area, the average annual precipitation is about 53.14 mm, while the annual potential evaporation is much higher, about 2196 mm. The average annual temperature is about 3.9 °C, which is typical of an inland arid climate. The basin covers an area of 1.02 <sup>×</sup> <sup>10</sup><sup>6</sup> km<sup>2</sup> and is the largest inland river basin in China

(Figure 1). In response to the global climate change and increasing human activities, the natural ecosystems in the basin are facing a series of crises and challenges. Its fragile ecological environment possesses abundant natural resources [24]. The drainage network in the basin consists of the main stream of the Tarim River and 144 drainage systems associated with nine major tributary basins: the Yarkand River, the Aksu River, the Kaidu-Kongque River, the Hotan River, the Kaxgar River, the Weigan River, the Dina River, the Keriya River, and the Qarqan River. The tributaries to the main stream of the Tarim River form a centripetal shape around the Tarim Basin [25]. The Tarim River is a dissipative inland river whose runoff is mainly supplied by meltwater from glaciers and snow. The sources of the Tarim River runoff include glacial melt water, accounting for 48.2%; precipitation in the form of rain and snow, accounting for 27.4%; and river base flow, accounting for 24.4% of the total [26,27]. The Tarim River Basin is a typical oasis agricultural production area in China; the oasis area in the basin has been increasing, and arable land has increased during the past 30 years. This has led to a shortage of water in the basin, which is mainly used for agriculture. The demand for agricultural irrigation water is large and accounts for about 96% of the total water use in the Tarim River Basin [28].

**Figure 1.** Map of the study area showing the Tarim River Basin and its nine main tributary river basins: the Kaidu-Kongque, Aksu, Weigan, Yarkand, Qarqan, Dina, Hotan, Kaxgar, and Keriya river basins; (**a**) the elevation of this area ranges from 781 to 8538 m above sea level; the spatial and temporal distribution of (**b**) temperature and (**c**) precipitation changes from 1990 to 2015. The black triangles represent the meteorological stations. Blue and red represent the increase and decrease, respectively, and the size of the triangle represents the magnitude of the change.

#### **3. Materials and Methods**

#### *3.1. Materials*

To analyze the dynamic evolution of the Tarim River Basin and its ecotone, and the factors controlling the observed changes, this article mainly uses remote sensing, land use, and meteorological data collected between 1990 and 2015. The analyzed historical trends in such parameters as the area of land-use transfer were subsequently used to extrapolate and predict the potential future changes in the study area.

#### 3.1.1. Remote Sensing Data

Landsat 5 Thematic Mapper (TM)/Landsat 7 Enhanced Thematic Mapper (ETM) +/Landsat 8 Operational Land Imager (OIL) satellite imagery from 1990, 2000, and 2015 (a total of 78 pieces, Table S1) was used in this study. The acquisition time was from June to September of each year. The cloud cover was less than 10%, and the pixel resolution of

the data set was 30 m × 30 m. The ENVI 5.3 software was used to perform radiometric calibration, atmospheric correction, and Normalized Difference Vegetation Index (NDVI) extraction calculations on the remote sensing images.

#### 3.1.2. LUCC Data

The land-use data set for the Tarim River Basin was obtained from the existing remote sensing monitoring data set of land use in China. It was provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (http://www. resdc.cn). In this paper, we used three periods of Chinese land-use data (30 m × 30 m), collected in 1990, 2000, and 2015. The database offers the most comprehensive coverage of China's land use/land cover and has been used in a number of published studies [29,30]. The land-use types were classified into six categories: arable land, forest, grassland, water, built-up land, and unused land (Table S2).

#### 3.1.3. Meteorological Data

Monthly scale meteorological information on temperature, precipitation, wind speed, humidity, and pressure from 1990 to 2015 for 26 meteorological stations in the Tarim River Basin were used to describe the recent changes in climatic conditions. The meteorological data were obtained from the China Meteorological Science Data Sharing Service Network, which have good continuity and has been tested for consistency. Station selection required the following: (1) the station was a national ground meteorological station, and (2) missing data accounted for less than 1% of the total data.

#### 3.1.4. Groundwater Data

We selected groundwater-level data from the Yarkand, Kaxgar, and Weigan river basins. Three groundwater monitoring wells were selected for each basin. Groundwater data for the Yarkand River Basin (2004–2010) were obtained from the Kaxgar Hydrological Bureau, and the groundwater data for the Weigan River Basin (2000–2012) and the Kaxgar River Basin (2004–2010) were obtained from groundwater monitoring wells deployed by the Xinjiang Institute of Ecology and Geography of the Chinese Academy of Sciences for multi-year actual measurements.

#### *3.2. Methods*

#### 3.2.1. NDVI Calculation

Vegetation indices are often used for vegetation analysis [31–33] and are typically formed by combining certain bands of image spectral data that possess vegetation-sensitive properties. Widely used vegetation indices include simple vegetation indices, ratio vegetation indices, NDVI, and transformed normalized vegetation indices [33,34]. The NDVI was used as the analysis index. It was calculated as follows:

$$\text{NIDVI} = (\text{NIR} - \text{R}) / (\text{NIR} + \text{R}) \tag{1}$$

where NIR is the TM near-infrared band value and R is the TM visible-red band value. The range of the image element values was -1 ≤ NDVI ≤ 1. Negative values indicate that the ground cover consists of clouds, water, snow, etc.; 0 indicates the occurrence of rock or bare soil; positive values indicate that there is vegetation cover, and the values increase with increasing coverage.

ENVI 5.3 was used to perform radiometric calibration, atmospheric correction, and mosaic and other processing of the Landsat images. The NDVI calculation tools were used to calculate and output images from the three selected years separately. The calculation results were then loaded into ArcGIS10.6, which was used to eliminate outliers. The GIS raster calculator was used to extract the desert-oasis ecotone.

Sun et al. (2020) indicated that when the NDVI is between 0.05 and 0.35, the scope of the delineated ecotone can match the actual scope of the ecotone to a greater extent [35]. On this basis, we also combined visual interpretation and field verification to exclude

the artificial protection forest and arable land in this range, and with this criterion, the desert-oasis ecotones of the Tarim River Basin in 1990, 2000, and 2015 were obtained.

#### 3.2.2. Land-Use Transfer Matrix

This research quantitatively studied the conversion between land-use types in different years. ArcGIS and MATLAB were used to process the land-use TIFF data from the study area and to calculate the land-use transfer matrix from 1990 to 2015.

The land-use transition matrix was derived from system analysis and was used to quantitatively describe the mutual feedback relationship between the system state and the state transition. Put differently, the matrix describes the change in conditions of the system from time T to time T + 1 [36]. At present, it is widely used to describe the internal characteristics of the transfer structure and transfer direction of land-use/landcover changes in basins. This method can not only describe the structural characteristics of the land-use area within a period of time but also effectively describes the transfer area and transfer direction of various land-use types at the beginning and the end of the period, as follows:

$$\mathbf{S}\_{\mathrm{ij}} = \begin{bmatrix} \mathbf{S}\_{11} & \cdots & \mathbf{S}\_{\mathrm{ln}} \\ \vdots & \ddots & \vdots \\ \mathbf{S}\_{\mathrm{n}1} & \cdots & \mathbf{S}\_{\mathrm{nn}} \end{bmatrix} \tag{2}$$

where S is the area; i and j (i, j = 1,2..., n) are the land-use type before and after transfer, respectively; *Sij* is the area where land use changes from type i to type j; and n is the number of land-use types before and after the transfer.

#### 3.2.3. The Standardized Precipitation Evapotranspiration Index

The Standardized Precipitation Evapotranspiration Index (SPEI) is widely used in global and regional drought detection and characterization [37–40]. This study used monthly site data and the Penman–Monteith formula to calculate differences in potential evapotranspiration, which were needed to calculate the SPEI. Specifics of the calculation can be found in Allen et al. [2].

A positive SPEI value indicates a relatively wet condition, whereas a negative SPEI value indicates a dry state. An SPEI value between -0.5 and 0.5 indicates normal condition. In this study, SPEIs at 3-month, 6-month, 9-month, and 12-month timescales were used for analysis. The ranges of SPEI values were divided into five classes according to the national meteorological drought scale (Table S3).

#### 3.2.4. The CA-Markov Model

The CA-Markov model was used herein to simulate and predict future land use in the Tarim River Basin. The Markov model is a method for predicting the probability of occurrence at a specified time based on the Markov chain process theory. It is often used for the prediction of geographic events with no after-effect characteristics [41,42]. The evolution of land use has the nature of the Markov process. In fact, land-use type corresponds to the "possible state" of the Markov process, and the area or ratio of the conversion between land-use types can be represented as a state transition probability matrix [43]. The model is expressed as follows:

$$S\_{(t+1)} = P\_{\vec{l}\vec{j}} \times S\_{\vec{l}} \tag{3}$$

$$P\_{ij} = \begin{bmatrix} P\_{11} & \cdots & P\_{1n} \\ \vdots & \ddots & \vdots \\ P\_{n1} & \cdots & P\_{nn} \end{bmatrix} \tag{4}$$

$$\{0 \le P\_{ij} < 1 \text{ and } \sum\_{j=1}^{n} P\_{ij} = 1 (i.j = 1, 2 \cdots n)\}\tag{5}$$

In Equation (3), S(t+1) is the system state at t + 1 and Pij is the state transition probability matrix.

The cellular automata (CA) model is a lattice dynamics model with discrete time and space states. It focuses on the interaction of different temporal and spatial characteristic cells and has powerful spatial calculation and simulation capabilities [44,45]. In terms of land-use prediction, the Markov model focuses on the prediction of the amount of land-use change but it cannot spatially express the spatial distribution of the various types of land-use changes. The cellular automata model can express the spatiotemporal dynamic evolution of complex spatial systems, thereby making up for the deficiency in the Markov model [46].

This study applied the CA-Markov model in the Idrisi 17.0 software to land-use data collected in 1990, 2000, and 2015. The first step in the approach used the 2000 data as the starting year and predicted the land use in 2015. The land-use predictions were then compared with the actual land use in 2015 to verify the reliability of the CA-Markov model simulation. Once verified, the spatial patterns in land use in 2030 were predicted using 2015 as the starting year (Figure 2). The images involved in the processing of the Idrisi software are portrayed as raster data, and the land raster size used in this analysis was 30 × 30 m. The spatial data processing was completed using ArcGIS software.

**Figure 2.** Flow chart showing primary steps in the simulation process.

#### 3.2.5. The Kappa Index

The Kappa index is often used to interpret remote sensing accuracy and to evaluate the similarity of two spatial maps. The Kappa index was used in this study to verify the accuracy of the CA-Markov model for simulating the evolution of land use in each tributary basin [47]. The index is calculated as follows:

$$kappa = \frac{P\_o - P\_c}{P\_p - P\_c} \tag{6}$$

$$P\_o = \frac{n\_1}{n}, P\_c = \frac{1}{N} \tag{7}$$

where *P<sup>o</sup>* is the proportion of the raster that is correctly simulated, *P<sup>c</sup>* is the desired proportion of correctly simulated raster grid cells, *P<sup>p</sup>* is the proportion of correctly simulated grids for an ideal classification, n is the total number of grids, *n*<sup>1</sup> is the number of grids that are correctly simulated, and N is the number of land-use types (N = 6 in this study). The degree of consistency is weak when the Kappa index is less than 0.2 and significant when the Kappa index is greater than 0.8. A higher Kappa index means a better model simulation. (The detailed relationship between the Kappa index and consistency is shown in Table S4.)

#### **4. Results**
