Remotely-Sensed Identification of a Transition for the Two Ecosystem States Along the Elevation Gradient: A Case Study of Xinjiang Tianshan Bogda World Heritage Site

Abstract

The alpine treeline, as an ecological transition zone between montane coniferous forests and alpine meadows (two ecosystem states), is a research hotspot of global ecology and climate change. Quantitative identification of its elevation range can efficiently capture the results of the interaction between climate change and vegetation. Digital extraction and extensive analysis in such a critical elevation range crucially depend on the ability of monitoring ecosystem variables and the suitability of the experimental model, which are often restricted by the weak intersection of disciplines and the spatial-temporal continuity of the data. In this study, the existence of two states was confirmed by frequency analysis and the Akaike information criterion (AIC) as well as the Bayesian information criterion (BIC) indices. The elevation range of a transition for the two ecosystem states on the northern slope of the Bogda was identified by the potential analysis. The results showed that the elevation range of co-occurrence for the two ecosystem states was 2690–2744 m. At the elevation of 2714 m, the high land surface temperature (LST) state started to exhibit more attraction than the low LST state. This elevation value was considered as a demarcation where abrupt shifts between the two states occurred with the increase of elevation. The identification results were validated by a field survey and unmanned aerial vehicle data. Progress has been made in the transition identification for the ecosystem states along the elevation gradient in mountainous areas by combining the remotely-sensed index with a potential analysis. This study also provided a reference for obtaining the elevation of the alpine tree line quickly and accurately.

1. Introduction

As the ecological transition zone of montane coniferous forests and alpine meadows, the alpine tree line is extremely sensitive to environmental factors such as temperature, precipitation, and soil, which has great ecological vulnerability [1]. With the change of environmental climatic factors, the vegetation distribution show a changing or narrowing position of the alpine tree line [2,3]. The Fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change (IPCC) noted that recent climate changes have generated widespread impacts on human and natural systems. The atmosphere and oceans are getting warmer, and the amounts of snow and ice are decreasing. Many terrestrial species have shifted their geographic ranges, abundances, and species interactions in response to ongoing climate change [4]. As a “monitor” of climate change, the alpine tree line is sensitive to climate change and can help humans quickly understand the interaction between climate change and vegetation [5,6,7]. In addition, the ecosystem is facing unprecedented pressure, which will restrict regional sustainable development in different regions. Consequently, it is one of the hotspots to study the influence of climate change by obtaining the elevation range of the alpine tree line. At the same time, strengthening the study of the alpine tree line is significant for the protection and sustainable development of Natural Heritage Sites.

Many researchers have carried out relevant studies on the alpine tree line related to the formation and restriction mechanisms, location fluctuations under climate change, vegetation community diversity patterns, and demarcation elevation identification. Some studies have taken the limitations of the alpine tree line formation as the entry point and used the experimental methods of physiology and ecology to conduct studies on temperature as the main controlling factors [2,8,9,10,11,12]. Precipitation, propagation limitation, carbon limitation, environmental stress, and external disturbances have been analyzed as the causes of formation for the alpine tree line [3,9,12,13,14,15,16,17]. The direct effects of climate change with stored terrestrial carbon include high temperatures, drought, and windstorms. The indirect effects include increasing risks of fires, pests, and disease outbreaks. It is projected that the rise in tree mortality will occur in many regions during the 21st century [4]. The relationship between alpine tree line distribution and climate change has been analyzed by means of tree chronology methods such as tree rings, sporopollenin, and sediment magnetic susceptibility, paleontology, and isotope studies [13,18,19]. The demarcation elevation and population structure are two important characteristics of the alpine tree line [20,21]. In addition, the alpine tree line was taken as an ecological transition, and the position changes of which were extracted. In the early stages, methods such as ecological field surveys [11,22,23] and historical data collection and collation [2,24] were usually used to obtain the demarcation elevation of the alpine tree line. However, due to the influence of environmental factors such as topography and climate, long cycles, and the high cost of human and material resources, there are many areas that cannot be accessed by humans in high mountains, which limits the feasibility and the objectivity of the research. The validation sites in situ have certain regional limitations. The historical data cannot accurately reflect long-term climate change and the current influence of human activities on the alpine tree line. With the development of spatial information technology such as 3S (RS, remote sensing, GNSS, global navigation satellite system, GIS, geographic information system), it is possible to observe a long-time series of the macroscopic area represented by RS technology in order to deal with multi-source data easily. It can provide a broad platform for studying the demarcation elevation of the alpine tree line. RS and GIS technologies have been introduced into the study of the alpine tree line. Remotely-sensed classification information and digital elevation models (DEM) were used to extract the demarcation elevation of the alpine tree line, which greatly improved the efficiency and comprehensiveness of the study of alpine tree line locations [21,25,26,27,28]. However, there are still some problems in the study of demarcation elevation. First, although the elevation location of the alpine tree line can be obtained based on remotely-sensed classification information and DEM, it relies too much on the classification information as well as the topographic factors and ignores the analysis of the ecosystem state. In fact, the alpine tree line is the area where the ecosystem state shifts. With uncertainties such as threshold window selection, the accuracy of the extraction results needs to be further improved [21,27]. Second, based on statistical methods and the stable-state principle, the data collected in situ were used to analyze the elevation distribution of larch forests, but this study did not have any relevant results about the elevation of alpine tree line, nor did it use RS and GIS technology [29]. Lastly, the highest, lowest, and mean elevations of the alpine tree line ecotone were extracted in the study [25]. However, there is a lack of field validation data, and the accuracy of results needs to be verified [25].

Since the 19th century, the alpine tree line has been often considered to be consistent with the 10 °C isotherm of the warmest month, but this conclusion applies only to the Alps and the Colorado Rocky Mountains [8,9]. The tree line of Mount William in New Guinea (3850 m, 5°S) is located in a region with an average temperature of 5.6 °C, while the tree line of the Craigieburn Mountains in New Zealand (1300 m, 43°S) is in line with the 11.6 °C isotherm of the warmest month [2]. In order to analyze the relationship between the alpine tree line and a universal temperature, Körner [2] discussed the causes of formation for the tree line from a global point of view, and noted that the location of the tree line was closely related to the average temperature during the growing season. Subsequently, Körner and Paulson [11] used the root temperature records of 46 tree line points ranging from 68°N to 42°S in 1996–2003, and found that the tree line in high altitude climates correlated with the average temperature of 6.7 °C during the growing season. For the average temperature of the growing season in different climatic zones, there was a narrow floating range of 2.2 K. It indicated that the growth of alpine tree line plants was controlled by temperature and there existed a low temperature threshold. Jobbágy and Jackson [10] utilized a global climate database and 115 points of alpine tree line data (excluding tropical zones), and the statistical analysis showed that the change in the tree line was correlated with an annual average temperature and a seasonal heat radiation value. There is not enough solar radiation data in most parts of the world, but, fortunately, the temperature of the canopy near the tree line is relatively close to the temperature of the air. Woodward and Körner [30,31] studied the leaf extension of different species at different elevations by means of an electronic assistant measuring instrument, and showed that the leaf extension of different species at different elevations was limited by temperature. Hence, the temperature is the main factor for controlling the formation of the alpine tree line, which has been supported by most studies. In the case of relatively moist conditions, the temperature is the main factor in the formation of the Bogda alpine tree line [32,33]. The temperature data is usually acquired from meteorological stations by point measurements. However, the distribution of meteorological stations in a mountainous area is extremely limited [34,35]. Generally, different interpolation methods are used to obtain continuous surface temperature data, and the results of different methods are also different, which greatly reduce the reliability of data. High-precision LST can be derived from remotely-sensed data. Therefore, we can make use of the Thermal Infrared Sensor (TIRS) data with a spatial resolution of 100 m to retrieve LST. LST, as one of the priority parameters of the International Geosphere-Biosphere Programme, can provide spatial and temporal information on the energy balance of the surface. It is key for remote sensing applications in drought prediction, crop water pressure estimation, plant physiological monitoring, crop yield estimation, numerical weather prediction, and climate prediction [36,37]. Thus, LST was selected as the index to extract the elevation range of the alpine tree line.

In view of these problems, we chose the Bogda component of the Xinjiang Tianshan World Natural Heritage Site (hereinafter referred to as Bogda) as our study area. By combining the 3S technology with a potential analysis model, the LST index was utilized to identify the transition for the two ecosystem states along the elevation gradient, and the results were validated. This study can provide some new ideas and methods for extraction and dynamic change-monitoring research of the alpine tree line.

2. Materials and Methods

2.1. Study Area

The Bogda is the eastern part of the Xinjiang Tianshan World Natural Heritage Site, located in Fukang County and Ürümqi City, in Xinjiang Uygur Autonomous Region (Figure 1). It is a wet island in the center of the arid desert belonging to the continental temperate climate zone. Most parts of the Bogda are sheltered by the windward slope intercepting currents from the Atlantic and Arctic Oceans. Because the mountainous body is a barrier, the temperate continental climate is not as clear in this area. It is warmer in the winter and cooler in the summer, with enough rainfall and relatively deep snow. The elevation gradually becomes higher. On the northern slope of the Bogda Peak, within a horizontal distance of 40 km, the elevation rises from 800 m to 5000 m, with six altitudinal natural zones. Montane coniferous forests and alpine meadows are distributed between 1600 m and 3300 m. The annual average temperature is 2.5 °C. In January, which is the month with the lowest temperatures, the average temperature is –12.4 °C. In July, which is the month with the highest temperatures, the monthly average temperature is 15.9 °C. The annual average precipitation is 444 mm [38].

The Bogda has a typical altitudinal natural zone, covering almost all the natural landscape types and ecosystems in central Asia. With the increase of the elevation, the altitudinal natural zones exhibited six types: temperate desert (800–1100 m), montane steppes (1100–1650 m), montane coniferous forests (1650–2700 m), alpine meadows (2700–3300 m), alpine cushion vegetation (3300–3700 m), and ice and snow (>3700 m) [38]. Within a horizontal distance of 40 km, the elevation rises from 800 m to 5000 m, with six altitudinal natural zones. It is rare that, within such a short horizontal distance, the integrated altitudinal zone is typical as well as representative of the temperate arid zone, which reinforces the global significance of the site in the study of the succession of montane ecosystems. Meanwhile, this area has also been accepted by the United Nations Educational, Scientific and Cultural Organization (UNESCO) as a Man and Biosphere Reserve [38]. In the study area, montane coniferous forests are dominated by Picea schrenkiana, while alpine meadows are composed of Kobresia capillifoli (Figure 2).

2.2. Data

2.2.1. Landsat Data

The Landsat-8 Level-1 data on July 28, 2016 with good vegetation cover was selected. The track number is 142/30 (The data set is provided by the Geospatial Data Cloud site, Computer Network Information Center, Chinese Academy of Sciences. http://www.gscloud.cn). On February 11, 2013, National Aeronautics and Space Administration (NASA) successfully launched the Landsat-8 satellite with two sensors, which include the OLI (Operational Land Imager) and the TIRS. The Landsat-8 Level-1 data products have undergone systematic radiation corrections and geometric corrections. The OLI instrument consists of eight multispectral bands (resolution 30 m) and one panchromatic band (resolution of 15 m). The resolutions of the thermal infrared bands 10 and 11 are 100 m, but the downloaded data in this study has been resampled to 30 m. Unlike previous Landsat series sensors, the TIRS instrument has two thermal infrared channels (bands 10 and 11) similar to those for bands 31 and 32 of the Moderate Resolution Imaging Spectroradiometer (MODIS) [39].

The Bogda is a typical mountain ecosystem. The radiation received by sensors is the result of complex interactions between the sun, the atmosphere, and the Earth’s surface. Therefore, in the mountainous areas, the radiation on different aspects has significant differences, which ensures the characteristics are reflected by remotely-sensed images different from the real features [40]. The influence of topography on the surface reflectance retrieval cannot be neglected [41]. This study performed topography corrections (SCS+C), atmospheric corrections, and cropping of Landsat OLI data (excluding panchromatic bands). The normalized difference vegetation index (NDVI) is a vegetation index, which is widely used to monitor vegetation growth and cover [42,43]. Then, NDVI was computed by the band 4(Red) and band 5(Near-Infrared, NIR) based on the processed image data, as shown in Formula (1).

$$NDVI=\frac{R_{NIR}-R_{Red}}{R_{NIR}+R_{Red}}$$

(1)

where $R_{NIR}$ and $R_{Red}$ are the spectral reflectances in the NIR and red bands, respectively.

2.2.2. Field Data

We developed field routes that were distributed on both sides of the Sangong River (the center of the results), and the manpower-accessible areas along the routes that were set up as survey points. From July 4, 2018 to July 14, 2018, unmanned aerial vehicles and high-precision Global Positioning System (GPS) devices were utilized to record the location information, elevation information, and a large number of field photos. Field routes and survey points are shown in Figure 3. There were 23 points in situ. Among them, seven points were used to validate the results, and the rest were used to record the land cover types along field routes.

2.2.3. Other Data

Other data were selected for analysis in the study (

Table 1. The information of other data.

Type	Resolution (m)	Sources
SRTM _1arc _v3	30	U.S. Geological Survey (USGS)
Sentinel-2 L1C	10	ESA data distribution website
Land cover data	30	Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC)

). The SRTM_1arc_v3 (SRTM, Shuttle Radar Topography Mission) are expressed in geographic coordinates (latitude/longitude) and are horizontally referenced to WGS-84 and vertically referenced to the EGM-96 Geoid, and meet the absolute horizontal and vertical accuracies of 20 m (circular error at 90% confidence) and 16 m (linear error at 90% confidence), respectively (https://lta.cr.usgs.gov/SRTM1Arc.data Description).

The Sentinel-2 Earth Observation Mission is part of European Space Agency’s (ESA) Copernicus Program, which consists of Sentinel-2A and Sentinel-2B with 13 bands. The spatial resolution of bands 2 (blue), 3 (green), 4 (red), and 8 (NIR) of the Sentinel-2 L1C data is 10 m, and the date of the used image was August 4, 2017. In addition, the land cover data for 2015 is generated by a manual visual interpretation based on the 2010 land cover data and the Landsat-8 remotely-sensed images.

2.3. Methods

2.3.1. Mono-Window Algorithm

The ground emissivity was derived from the OLI data, and the LST was retrieved from the TIRS 10 data. The approaches for LST retrieval mainly include the split-window algorithm, the single-channel method, and the multi-channel method. Landsat-5 Thematic Mapper (TM) and Landsat-7 Enhanced Thematic Mapper Plus (ETM+), with a single thermal infrared band, typically use a single-channel method to retrieve the LST. Landsat TIRS’s two thermal bands are designed for applications such as atmospheric correction and the retrieval of LST using split-window algorithms. However, the USGS had announced that the calibration parameters of band 11 are uncertain, and the split-window algorithm is not encouraged [44,45,46]. In addition, the accuracy of LST retrieved from band 10 was higher than that from band 11 [47,48].

Single-channel algorithms for the retrieved LST usually include the radiative transfer equation-based method, the mono-window algorithm developed by Qin et al. [49] and the generalized single-channel algorithm developed by Jiménez-Muñoz and Sobrino [50]. The mono-window algorithm has been applied to LST retrieval of Landsat series data [46,51,52,53]. Qin et al. [51] utilized the mono-window algorithm to retrieve the LST across the Israel-Egypt border with an accuracy of 1.1 °C. Hu et al. [52] obtained the LST in a suburb of China based on the mono-window algorithm and the generalized single-channel algorithm, and then validated the results through the measured data in situ. The results showed that the overall average errors were 0.83 °C and 1.08 °C, respectively. Thus, the mono-window algorithm is used to retrieve the LST.

The approach devised herein was based on the joint use of: (1) the spectral radiance that is calculated from the thermal radiance, (2) the brightness temperature at the satellite level can be calculated, (3) the ground emissivity is estimated according to the method [54,55,56], and atmospheric profile parameters are obtained from the atmospheric parameter calculation website provided by NASA, (4) blackbody brightness temperature is obtained based on previous parameters, (5) and the LST is retrieved. Furthermore, the main process is as follows [51].

The spectral radiance received by the sensor has the following relationship with its DN (Digital Number, pixel values of remotely-sensed image) value.

$$L_{\left(\lambda \right)}=L_{min\left(\lambda \right)}+\left(L_{max\left(\lambda \right)}-L_{min\left(\lambda \right)}\right)Q_{\left(dn\right)}/Q_{max}$$

(2)

where $L_{\left(\lambda \right)}$ is the spectral radiance received by the sensor and $Q_{\left(dn\right)}$ is the gray level. The DN value is converted to radiance by radiometric calibration.

The brightness temperature at the satellite level can be computed by using the following equation.

$$T_{10}=K_2/\ln (1+K_1/L_{\left(\lambda \right)})$$

(3)

where $T_{10}$ is the satellite brightness temperature of TIRS 10 and $K_1$ and $K_2$ are prelaunch calibration constants. For Landsat 8, $K_1=774.89W/\left(m^2·sr·um\right)$, and $K_2=1321.08$ K.

The LST is retrieved from the Landsat 8 data as follows.

$$T_s=\left[a\left(1-C-D\right)+\left(b\left(1-C-D\right)+C+D\right)T_{10}-D{T}_a\right]/C.$$

(4)

$$C=\epsilon ·\tau$$

(5)

$$D=\left(1-\tau \right)·\left[1+\left(1-\epsilon \right)·\tau \right]$$

(6)

where $T_s$ is the LST, $T_a$ refers to the atmospheric temperature, $\epsilon$ is ground emissivity, and $\tau$ is atmospheric transmittance. In the formula, $a$ and $b$ are the coefficients (a = −67.355351, b = 0.458606).

2.3.2. Total Shortwave Broadband Albedo

Surface albedo is defined as the ratio of the solar radiation reflected from Earth’s surface to the solar radiation incident upon it [57]. Land surface broadband albedo is a key biophysical variable in many scientific applications [58] (such as the surface energy balance [59] and global meteorological changes [60]), and it will change with the variations in vegetation type [61,62]. Higher-resolution albedo products derived from Landsat data can describe the heterogeneous environments and land cover changes at a small scale [63,64]. Liang [58] proposed an algorithm to estimate the land surface broadband albedo based on Landsat data under various atmospheric and surface conditions. This estimation approach can be used to derive the surface shortwave albedo from top-of-atmosphere (TOA) reflectance using a statistical relationship established by radiative transfer simulations. Then, ground measurement data of several cover types were used to validate the simulation results, and the average residual standard errors of the estimated broadband albedo were around 0.02 [65]. The linear equation is given as follows.

$${\alpha }_{\lambda }={\displaystyle \sum }{\rho }_i^{TOA}·c_i+c_0$$

(7)

where ${\alpha }_{\lambda }$ is the surface albedo for the spectral range of $\mathsf{\lambda }$, ${\rho }_i^{TOA}$ is the TOA reflectance for band i, and $c_i$ and $c_0$ are the regression coefficients [58,66]. The values of $c_i$ and $c_0$ are given in

Table 2. The values of the $c_i$ and $c_0$.

Band i	2	4	5	6	7
$c_i$	0.356	0.130	0.373	0.085	0.072	--
$c_0$	--	--	--	--	--	−0.0018

.

In this study, Liang’s algorithm was used to estimate the total shortwave broadband albedo (SBA) of the interest area, and then the potential energy of the SBA was analyzed.

2.3.3. Statistical and Frequency Distribution Analysis

For mountainous ecosystems, LST is affected by topographic factors such as elevation, aspect, type of land cover, soil moisture, incident radiation, and slope [67]. The influence of the aspect on the Tianshan Mountains in Xinjiang is remarkable. The temperature levels in different regions are very different [68], and the vegetation distribution on the southern and northern slopes is clearly different. The northern slope of the Bogda is more humid than the southern slope. Thus, different vegetation types grow on the southern and northern slopes at the same elevation, among which the montane coniferous forests zone is the most clear. With the change of elevation, the coniferous forests on the north-facing slope and the steppes on the south-facing slope are often alternately distributed. Therefore, our study took into account the influence of land cover types and aspects on the differences of LST distribution. The relationship between the aspects and the LST was analyzed by combination with the actual situation of the Bogda.

If a system has only one state, the frequency distribution usually forms a single-peak mode under stochastic disturbances and environmental variations. If there are other states, it will form a multi-peak mode [69,70]. A multimodal frequency distribution of states can result from environmental drivers or from the existence of the multiple states in the system [69,70]. These frequency distributions can reflect the approximate shape of the modal distribution around the transition of two states. Before establishing the probability density distribution, the LST values were arcsine transformed to approximate a normal distribution.

The AIC was proposed and developed by Hirotugu Akaike in their research information theory. It is a standard for estimating the goodness-of-fit of statistical models [71]. When there are several competing models, the model with the smallest AIC value given by the maximum likelihood estimate of the parameters is desired [72].

$$AIC=2k-2ln\left(L\right)$$

(8)

where k is the number of models and L is the likelihood function.

The BIC was proposed by Schwarz in 1978. Similar to the AIC, it is also used for model selection. The BIC value of each model is calculated, and the model with the smallest BIC value is desirable [73]. For overfitting, both the AIC and the BIC have introduced penalty terms related to the number of model parameters, and the penalty for BIC is greater than that for AIC. When the number of samples is excessive, it can effectively prevent the model complexity from being too large, which is caused by excessive model accuracy.

$$BIC=kln\left(n\right)-2ln\left(L\right)$$

(9)

where k is the number of models, $n$ is the number of samples, and L is the likelihood function. The $kln\left(n\right)$ penalty term can effectively avoid a dimensional catastrophe when the dimension is sufficiently large and the training sample data is relatively small.

The frequency distribution of LST within the interest area was analyzed by establishing an empirical histogram. Then, the expectation-maximization algorithm was utilized to find the normal distribution of the optimal matching number. The number of optimal modes can be determined through AIC and BIC, which is related to the number of states. The variation of states implies the potential of critical transitions [74].

2.3.4. Potential Analysis

As a state parameter, potential energy is also a kind of energy that is stored in a system. From the ecological point of view, the ecosystem should maintain the stability of existing structures and functions at a specific time and space scale. However, the study has revealed that such patterns, discovered through spatial gradients, can indicate different dynamic states [75]. In addition, ecosystems can be regarded as composed of multiple state modes, and each mode tends to present the most similar one with the adjacent state. Thus, when different state variables vary abruptly and exhibit an area of coexistence along the environmental gradient, a critical transition may occur between these states [29,76]. Potential analysis is a technique, which is used to establish a link between environmental gradients and state variables. Therefore, the potential analysis model was used to analyze whether the LST index of landscape states exhibited this pattern along the elevation gradient. The altitudinal vegetation differentiation can be treated as a transition process along the elevation gradient, which was caused by altitudinal differences in climate. Considering the transition of the ecosystem states along the elevation gradient as a nonlinear dynamic system, different landscapes (state variables) will shift along the elevation gradient (driver). The potential analysis can obtain potential energy values under different state variables by using a potential function, which consists of state variables and driving factors [69,75,77,78,79]. Thus, this method can calculate potential energy values of a landscape state along a driving factor gradient. Therefore, the elevation range in which two states co-occurred was a critical transition from one state to another in this study.

In a potential analysis, it is assumed that there exists a stochastic system with a potential function [69].

$$dz=-U\text{'}\left(z\right)dt+\sigma dw$$

(10)

where $U\left(z\right)$ is the potential function, $z$ is the state variable (in this case, LST), $\sigma$ is the noise level, and $\mathrm{d}w$ is a noise term. Given the relationship between the potential energy and the stationary probability density, the Fokker-Planck equation is applied to link them [75]. The formula for the potential energy of the system is as follows.

$$U=-\frac{{\sigma }^2}{2}log\left(Pd\right)$$

(11)

where $Pd$ is the probability density function of the state variables (in this case, LST). The probability density of LST was estimated using the ksdensity function in MATLAB. The bandwidth that controls the smoothness of the estimator is calculated as 1.06s/n^1/5 (s and n are the standard deviation and the length of the data set, respectively; see References [69,76,77]).

The abrupt nature of a shift in LST states along the elevation gradient was studied by tracking the local minimum potential energy (local minimum potential energy correspond to the different states [69,76,77,79]). In the ascending order of elevation, using the moving window (the window size was 100 [76]) to calculate the potential energy and local minimum for each iteration. Then, the potential energy and the number of states shift with increasing elevation were depicted by mapping the variation of LST along elevation gradients. We attempted to figure out whether the LST varies continuously or abruptly with an increasing elevation gradient. If the LST had a sudden change, the critical transition of the two states occurred at an elevation range. We defined this range as the demarcation elevations.

3. Results and Analysis

3.1. Relationship Analysis for the LST and Aspect

The shifts in land cover types with elevation (altitudinal natural zones) is the outstanding universal value (OUV) of the Bogda [38], which is affected by the spatial variability of water resources and heat energy in elevation and aspect. The Bogda is located in the eastern part of the northern Tianshan Mountains, and the northern Tianshan Mountains are in an east-west zonal mountainous system. In this study, except for the aspect analysis, we divided the area of interest into south-facing slope and north-facing slope. In addition, in aspect analysis, the aspect was classified into eight categories: the north (0° ~ 22.5°, 337.5° ~ 360°), northeast (22.5° ~ 67.5°), east (67.5° ~ 112.5°), southeast (112.5° ~ 157.5°), south (157.5° ~ 202.5°), southwest (202.5° ~ 247.5°), west (247.5° ~ 292.5°), and northwest (247.5° ~ 337.5°).

The vegetation types that grow on the southern-facing and northern-facing slope of the Bogda have significant differences. The montane coniferous forests mainly grow on the northern-facing slope, and the aspect distribution map is illustrated in Figure 4a (extracted from SRTM DEM data). Furthermore, the LST was retrieved by the mono-window algorithm (Figure 4b). Based on land cover types (Figure 2a), the minimums, maximums, means, and standard deviations of the LST for different land cover types were compared (Figure 4c). In the summer, the area with the highest LST (mean approximately 40 °C) was distributed in the low-covered and medium-covered grasslands corresponding to the temperate desert and the desert grassland, respectively. The area with lower LST (mean approximately −12 °C) was distributed in the snow and glaciers.

The woodland in Figure 2a was montane coniferous forests, the minimum and maximum of the LST were 15 °C and 34 °C, respectively, and the standard deviation was approximately 2.87. This is because the transpiration of the forest makes the temperature fluctuation relatively small. The high-covered grassland in Figure 2a was montane steppes and alpine meadows. The minimum and maximum of the LST were 3 °C and 40 °C, respectively, and the standard deviation reached approximately 3.85 °C. The main reason for the larger standard deviation was that the grassland has much wider distribution of elevation (1100–3300 m), and the elevation rise has a clear effect on the LST.

In order to study the relationship between the LST and aspect, the maximums, minimums, and means of the LST in eight aspects were calculated. The results, as illustrated in Figure 5, showed an insignificant fluctuation of the omnidirectional LST maximums, which, in the aspects of east, southeast, south, and southwest were slightly higher. Relatively, the fluctuation of minimums was more clear in the aspects of west, northwest, north, and northeast. The LST mean gradually increased from the east to the southeast, and then gradually decreased to the lowest on the northwest. In one word, the whole process exhibited some characteristics of the "sinusoidal" curve. The highest difference between the LST means in eight aspects can reach 6.43 °C.

The solar radiation and regional airflow can lead to different LST features. On one hand, the south-facing slope can receive more solar radiation, which shows higher LST. On the other hand, the regional airflow can affect the LST through the generation of land cover type. Influenced by Atlantic and Arctic currents, the climate in north-facing slope was wetter than in the south-facing slope, which makes montane coniferous forests the major land cover type. The larger heat capacity of the montane coniferous forests had a smaller standard deviation than that of montane steppes. The overpass time of the satellite was at local time 12:49 corresponding to 04:49 in GTM. At that time, the steppes heated up faster than the coniferous forests, so that the LST of the steppes also had a higher value.

Based on Landsat OLI data, the NDVI values were estimated (Figure 6). The high-value areas of NDVI were mainly distributed with montane steppes, montane coniferous forests, and alpine meadows, and the corresponding elevation was approximately 1600–3300 m, which was also the area of interest in this study. The NDVI of temperate desert and alpine cushion vegetation were relatively lower, and the NDVI of glaciers and snow were negative.

In order to extract the critical transformation along the elevation gradient between montane coniferous forests and alpine meadows, the influence of aspect on land cover types was considered. Within the elevation of 1600–3300 m, the land cover types of the north-facing slope shifted from meadow steppes to montane coniferous forests and then to alpine meadows while the land cover type of south-facing slope had no clear change. This distribution had been validated by field photos (Figure 7). Thus, the demarcation elevations of transition from the montane coniferous forests to the alpine meadows under the premise of eliminating the influence of the south-facing slope were researched and analyzed.

3.2. Critical Transitions in Altitudinal Zonality

3.2.1. LST Revealed the Presence of Two States

The results (Figure 8a) showed that the frequency distribution of LST was strikingly bimodal (at LST of approximately 20 °C and 28 °C), which indicates that there were distinct alternative modes. At the same time, the AIC and BIC of one to five normal distribution modes were simulated (Figure 8b). When the number of system modes changed from 1 to 2, both the AIC and BIC decreased rapidly. As the number of system modes increased, the AIC fluctuated slightly but remained stable and the minimum of BIC lay at two modes. Both indicated that LST had two modes in the study area. However, if the ecosystem has multiple states, there is a possibility of critical transition from one state to the other state [80,81]. Montane coniferous forests and alpine meadows were the only two possible states in the area of interest in July.

3.2.2. Detecting the Elevation Range of the Transition

The critical transition of two stable states along altitudinal zonality was calculated with elevation as the driver by using potential analysis techniques. The results of potential analysis showed details in how the two LST modes changed along the elevation gradient (Figure 9). The two models of LST were represented by two main land cover types in the area of interest. Moreover, one state corresponded to low LST (montane coniferous forests) and the other state corresponded to high LST (alpine meadows).

Below the elevation of 2690 m, there was a single state at low LST values (Figure 9b). Two states started to co-occur persistently at the elevation range from 2690 m to 2744 m (Figure 9c–e). Above the elevation of 2744 m, there was another single state at high LST values again (Figure 9f). Clearly, within the elevation of 2690–2714 m where two LST states co-occurred, the low LST state exhibited more attraction than the high LST state (Figure 9c). However, the low LST state started to exhibit less attraction than the high LST at the elevation of around 2714 m (Figure 9d). This trend continued until the elevation of 2744 m (Figure 9e). A sharp, discontinuous abrupt increase was observed around the elevation of 2744 m. In this process, the elevations (i.e., 2690 m, 2714 m, and 2744 m) of the transition for two ecosystem states were obtained.

3.3. Field Validation

The results showed that the transition’s elevation range was 2690–2744 m, and the demarcation elevation of two states shifting was 2714 m. It was consistent with the research conclusions of the Tianshan altitudinal zones given by other studies [38,82,83]. However, there were probable errors in estimating the key parameters, and this may be unavoidable. Therefore, field validation was conducted at the Bogda on 4–14 July 2018. Along the two sides of the Sangong River, the location information, elevations, and a large number of photos in situ were recorded using unmanned aerial vehicles and high-precision GPS devices (Figure 10b–d). In addition, Sentinel-2 L1C and Google Earth data were used to aid in validating the results (Figure 10a,e).

The validation of the results was provided by combining with the validation sites in situ, as shown in

Table 3. The validation of the results.

	Transition Range (Elevation/m)				State Shifting (Elevation/m)
	PtID	Starting (m)	Ending (m)	Difference (m)	Demarcation (m)	Difference (m)
Identified results		2690	2744		2714
Validation data	1	2664		−26
	2	2662		−28
	3				2718	+4
	4				2705	−9
	5		2732	−12
	6		2745	+1
	7		2757	+13

Note: “-” denotes when the validation sites in situ were lower than the identified results; “+” denotes when the validation sites in situ were higher than the identified results.

. The demarcation elevation differences of the critical transition were +4 m and −9 m, respectively (correspond to PtID 3–4), and the identified results were consistent with the validation sites. At the same time, the differences of the starting elevation were −26 m and −28 m, respectively (corresponding to PtID 1–2). However, the elevation of the field validation sites was lower than that of the identification results, which may be related to the location of the validation sites. PtID 1–2 were close to the Sigong River Basin, and the differences were slightly larger due to the influence of microtopography. In addition, the differences of the ending elevations were12 m, +1 m, and +13 m (corresponding to PtID 5–7), and the situation was almost the same as those of demarcation elevations. Consequently, it was shown that the extraction results were generally acceptable.

3.4. Potential Analysis of the SBA

The SBA spatial distribution was shown in Figure 11. The higher SBA values were mainly distributed with glaciers and snow (above the elevation of 3300 m). In the elevation range of 1600 to 3300 m, the SBA’s value of montane coniferous forests was significantly lower than that of meadow steppes. The value of SBA was relatively high in temperate desert regions (below the elevation of 1600 m). The same area of interest as LST covered was selected to analyze the distribution of the potential energy value of the SBA along the elevation gradient (Figure 12).

The results of potential energy for the SBA along the elevation gradient were shown in Figure 12. From Figure 12, it was found that there were two states of high albedo values and low albedo values. Below the elevation of 2740 m, where only the low albedo state was observed, most of the albedo was less than 0.1. Above the elevation of 2740 m, no sharp and discontinuous abrupt shift occurred in the low albedo state, but there was a significant fluctuation with a rapid rise to around 0.15, closed to the high albedo state. This elevation coincided with the ending of the transition range of the LST. During the elevation range when the two albedo states co-occurred (2740 m–2870 m), the low albedo in this scenario was significantly higher than the previous low albedo (below 2740 m), but slightly lower than the high albedo.

4. Discussion

It has been reported that albedo can be affected not only by fractional vegetation cover, but also by structural properties of vegetation such as the crown size, height, and shape of trees [84,85]. Combined with the actual situation, as shown in Figure 13, the area between the starting elevation and the ending elevation was the transition zone from coniferous forests to alpine meadows. Although coniferous forests and meadows co-occurred in this area, the number of coniferous forests was gradually decreasing. At the same time, the albedo increased slowly from below 0.1 to 0.1, but generally turned to the low albedo state of coniferous forests. However, above the ending elevation, there were mainly alpine meadows and a few shrubs (Figure 14). In addition, the albedo of the shrubs is higher than that of the alpine meadows but lower than that of the coniferous forests [76,86], and this may explain why the albedo states can show such a trend in Figure 12.

Therefore, elevation is a key factor in the altitudinal structure of vegetation in mountainous ecosystems [29]. Regardless of the metric considered, the presence of two states along the elevation gradient have been found. Through field validation, it is clear that the LST index can better reveal the transition from a mountainous coniferous forest to an alpine meadow. Many researchers have confirmed that temperature is the main factor controlling the formation of the alpine tree line by using discrete survey data or experimental data [2,9,10,12]. Our results were consistent with these studies.

Traditional research on the alpine tree line has mostly focused on field survey data or literature collection data [2,11,22,23,24]. However, on the spatial scale, if these discrete data are used to represent large mountainous systems, it will produce negative effects on the accuracy of mathematical simulation. On the temporal scale, it is difficult to form long-term data sets by a ground survey in such steep terrain because it would consume excessive human and material resources. In the past, there were few studies identifying ecosystem transitions along the elevation gradient by combining the remotely-sensed technology with the potential analysis [21,25,26]. For example, in the Bogda, Ji et al. [87] used the threshold analysis of the quantity ratios to obtain the demarcation elevation between montane coniferous forests and alpine meadows based on a scatterplot, which integrated the land cover classification (Landsat OLI data, on 28 July 2016) with the NDVI and DEM. The comparative analysis of the two studies was shown in

Table 4. Comparative analysis of the results.

	Elevation (m)
	Demarcation	Starting	Ending	Difference (Demarcation) PtID 3–4
Ji’ results [81]	2730	--	--	−12, −25
Our results	2714	2690	2744	+4, −9

. Comparing the two studies, we used the potential energy model to quantitatively identify the starting and ending elevations of the transition from the mountain coniferous forests to alpine meadows and the demarcation elevation of the shift for the two states. However, in Ji’s study, the results by the probabilistic statistical analysis depend on the size of the rolling window and threshold. It is impossible to generate a transition range using their methods. Moreover, our results were better than those of Ji, according to the validation in situ.

The tree line detection is more complicated than the direct division of the elevation range (e.g., aspect, slope, substrate, facilitation mechanism, microtopography, microclimate condition, local wind, distance to ridges, seed source, and disturbances). Based on the transition of ecosystem states, our method has made little progress in the identification of the alpine tree line location, and it is necessary to deeply explore the facilitation mechanism between different ecosystem states and the disturbances caused by the climate and by humans. In some areas, the alpine tree line has been shifted by anthropogenic impact [88] and global warming [17,89]. We have been to the Bogda for data acquisition and materials collection in 2017 and 2018, respectively. However, at present, we aimed to apply the method to identify the transition for the alpine tree line, without involving the change of the alpine tree line from the past to the present in this manuscript. Since remote sensing technology has gradually become an effective means to study a critical transition [78,90], exploring remotely-sensed indicators that have clear mechanistic links to the mountain ecosystem will facilitate the study of the altitudinal natural zone in the future.

The Bogda is a typical representative of the mountainous ecosystem in temperate arid zones and has typical altitudinal natural zones in these areas [91]. The distribution of altitudinal natural zones reflects the relationships between the biodiversity of lofty mountains in temperate arid zones and the changes in elevation, aspect, and slope. Therefore, it becomes a sensitive indicator of climate change in arid regions [38]. Based on this research, we aim to catch a breakthrough point for identifying the location distribution of the altitudinal natural zone in the Tianshan Mountains of Central Asia, and explore the relationship with global climate change. Furthermore, this study can provide a scientific basis for the sustainable development of the World Natural Heritage Site.

5. Conclusions

In this study, by integrating the remotely-sensed indicator with frequency distributions, AIC/BIC, and potential analysis, we identified two ecosystem states in the area of interest. The results indicated that the potential energy analysis model can be utilized to identify the number of states in the nonlinear dynamic system along the elevation gradient. The two ecosystem states co-occurred within the elevation range of 2690–2744 m, which is a transition from montane coniferous forests to alpine meadows. The shifts of the ecosystem states along the elevation gradient started at the elevation of 2714 m. There was also an abrupt and discontinuous critical transition at the elevation of 2744 m. The ecosystem states identification and potential analysis can be applied to the study of the altitudinal zonality in alpine areas. By detecting the critical transitions, the elevation range of each altitudinal zone can be identified.

Through the two indicators of LST and SBA, the shifts of different ecosystem states were revealed, but the transition details of separate indicators were different. Our results demonstrate the applicability of remotely-sensed indicators for identifying mountain ecosystem states, especially those that have a mechanistic link to ecosystem states, for exploring a transition of the states along the elevation gradient and improving the ability to study the location changes of the alpine tree line.