Dynamics and Drivers of the Alpine Timberline on Gongga Mountain of Tibetan Plateau-Adopted from the Otsu Method on Google Earth Engine

: The alpine timberline, an ecosystem ecotone, indicates climatic change and is tending to shift toward higher altitudes because of an increase in global warming. However, spatiotemporal variations of the alpine timberline are not consistent on a global scale. The abundant and highest alpine timberline, located on the Tibetan Plateau, is less subject to human activity and disturbance. Although many studies have investigated the alpine timberline on the Tibetan Plateau, large-scale monitoring of spatial-temporal dynamics and driving mechanisms of the alpine timberline remain uncertain and inaccurate. Hence, the Gongga Mountain on the southeastern Tibetan Plateau was chosen as the study area because of the most complete natural altitudinal zonation. We used the Otsu method on Google Earth Engine to extract the alpine timberline from 1987–2019 based on the normalized di ﬀ erence vegetation index (NDVI). Then, the alpine timberline spatiotemporal patterns and the e ﬀ ect of topography on alpine timberline distribution were explored. Four hillsides on the western Gongga Mountain were selected to examine the hillside di ﬀ erences and drivers of the alpine timberline based on principal component analysis (PCA) and multiple linear regression (MLR). The results indicated that the elevation range of alpine timberline was 3203–4889 m, and the vegetation coverage increased signiﬁcantly ( p < 0.01) near the alpine timberline ecotone on Gongga Mountain. Moreover, there was spatial heterogeneity in dynamics of alpine timberline, and some regions showed no regular trend in variations. The spatial pattern of the alpine timberline was generally high in the west, low in the east, and primarily distributed on 15–55 ◦ slopes. Besides, the drivers of the alpine timberline have the hillside di ﬀ erences, and the sunny and shady slopes possessed di ﬀ erent driving factors. Thus, our results highlight the e ﬀ ects of topography and climate on the alpine timberline on di ﬀ erent hillsides. These ﬁndings could provide a better approach to study the dynamics and formation of alpine timberlines.


Introduction
Structure, ecosystem processes, and functions of global vegetation communities and human livelihood have been seriously impacted adversely by climate change over the past 100 years [1]. Terrestrial ecosystems are undergoing long lasting, complex, and unpredictable changes as a consequence of variations in the global climate [2,3]. Vegetation, as a natural bond between soil, air and water, is crucial for global biodiversity conservation, ecosystem services and vision of human beings, however, it is highly sensitive to climate change [4]. Over the past few decades, terrestrial vegetation activity and growing season time in the Northern Hemisphere have increased with increase in temperatures [5][6][7]. Several cases worldwide have indicated that diverse vegetation patterns are greatly altered by global warming. For instance, the Arctic tundra has been encroached upon by shrubs, affecting ecosystem function, carbon storage, wildlife habitats, and hydrological status and processes [8,9]. Particularly, in the ecosystem ecotone (such as alpine timberlines), the vegetation community structure, function and spatial pattern are especially sensitive to global change [10,11].
The "alpine timberline" is a transition zone alternating between trees, shrubs, and grasslands along altitude gradients [12,13]. Environments in alpine timberline ecotone are characterized by lower annual temperatures and insufficient precipitation, resulting in which the trees experience cold climate and physiological drought that hamper their growth [14]. Additionally, the alpine timberline fluctuates with variations in multiple factors, such as dominant tree species, topography, and climatic conditions [15]. Alpine timberline dynamics in different regions currently present distinct characteristics. For example, based on nearly 70 years of historical maps and modern satellite imagery in the Ukrainian Carpathians, researchers found that the altitude of timberline had risen by 1.14 m year −1 , and the area covered by one-third of the pastures had increased [16]. In the alpine region of the Himalayas, the timberline shifted 10 m year −1 from the 1970s to 2006 [17]. Near the Yamunotri watershed, the alpine timberline moved to a higher altitude at an average rate of 4 m year −1 during the past 20 years [18]. However, the response trend of alpine timberline displacement to climate warming remains poorly understood, especially in the brittle ecosystem of the Tibetan Plateau.
The alpine timberline-critical for maintaining alpine biodiversity and functioning montane ecosystems-is the long-term result of natural environmental conditions and is extremely sensitive to global change [3,19]. The gradual decrease in temperature with elevation in montane systems greatly influences vegetation metabolic processes and alpine timberline formation [3,20]. Therefore, temperature is generally regarded as the primary driver of alpine timberline [21]. On a global scale, the alpine timberline possessed a trend consistent with the global ground isotherm (6.7 ± 0.8 • C) [22]. However, changing temperatures in different regions create an asynchronism between the dynamics of alpine timberline and global warming [23,24]. Studies in the European Alps indicate that alpine timberlines are strongly correlated with temperatures in the hottest month [13]. In northern Russia, summer temperatures and the vegetation growing season largely determine the location of the alpine timberline [25]. In China, the temperature during the growing season is the primary climatic factor affecting the altitude of the alpine timberline [26]. Additionally, rainfall is also an important factor regulating alpine timberlines in some regions and can improve growth conditions [27,28]. In short, water and heat conditions are the key drivers affecting alpine timberlines; however, the dynamic mechanisms involved remain undetermined as yet.
Gongga Mountain, located at the south-eastern edge of the Tibetan Plateau, belongs to the transition zone between the Sichuan Basin and the Tibetan Plateau. Because of the intact altitudinal zonation, Gongga Mountain is an ideal place to address the spatiotemporal dynamic and formation mechanism of the alpine timberline [29,30]. Previous studies obtained alpine timberline information mainly from historical data and field investigation [31,32]. Although field survey has benign accuracy and reliability, the relatively smaller scale restricts operational efficiency and broad applicability. Additionally, field surveys are inconvenient for monitoring long-term dynamic changes of alpine timberline on a large scale. Thus, exploring a new method to resolve dynamic tendency, spatial patterns, and drivers of alpine timberline is essential.
Remote sensing and geographic information technologies have recently become a practical means to identify the spatial distribution of alpine timberline [33,34]. The extraction methods of alpine timberline are mainly based on the supervised classification extraction and empirical normalized difference vegetation index (NDVI) values [35][36][37][38]. However, supervised classification requires the selection of training samples and human intervention, and these are characterized by several uncertainties with regard to the long-term monitoring of alpine timberlines. Meanwhile, the alpine timberline obtained by vegetation type boundaries is inaccurate, violating the timberline ecotone definition. Furthermore, empirical NDVI value is not very adaptable to different times, regions, and image data. The Otsu algorithm (the adaptive threshold algorithm) is introduced into Google Earth Engine to automatically determine the NDVI threshold of the alpine timberline ecotone, which addresses the aforementioned complications and improve computing efficiency. The Otsu algorithm finds an optimal threshold to distinguish relatively homogenous objects by maximising the inter-class variance [39][40][41]. Since the alpine timberline is an ecotone, the ecological vulnerability here is the lowest and the vegetation changes are the most unstable. Therefore, we retain vegetation and remove other land types (such as snow mountains, bare land, lakes and so on), and use the Otsu algorithm to obtain the NDVI threshold in the alpine timberline and to determine the location of the alpine timberline. The Google Earth Engine platform provides common remote sensing image data sets that can be called and calculated on online, providing convenience for the long-term monitoring of alpine timberlines [42,43].
Given these problems and new method, we calculated the NDVI based on Landsat imagery from 1987-2019, applied the Otsu method to Google Earth Engine, and obtained climate raster data by spatial interpolation that was based on terrain and generated topographic data by spatial analysis in the ArcGIS 10.2 software (ESRI, Inc., Redlands, CA, USA) across Gongga Mountain. Our primary objectives were: (1) determining the NDVI threshold of alpine timberline via the Otsu algorithm, (2) analysing spatial and temporal dynamics of the alpine timberline and exploring the effects of topography on alpine timberline distribution, and (3) revealing the comprehensive relationships among alpine timberline, topographic and climatic factors and discussing the variations of alpine timberline in future.
Remote Sens. 2020, 12, x FOR PEER REVIEW  3 of 19 alpine timberline obtained by vegetation type boundaries is inaccurate, violating the timberline ecotone definition. Furthermore, empirical NDVI value is not very adaptable to different times, regions, and image data. The Otsu algorithm (the adaptive threshold algorithm) is introduced into Google Earth Engine to automatically determine the NDVI threshold of the alpine timberline ecotone, which addresses the aforementioned complications and improve computing efficiency. The Otsu algorithm finds an optimal threshold to distinguish relatively homogenous objects by maximising the inter-class variance [39][40][41]. Since the alpine timberline is an ecotone, the ecological vulnerability here is the lowest and the vegetation changes are the most unstable. Therefore, we retain vegetation and remove other land types (such as snow mountains, bare land, lakes and so on), and use the Otsu algorithm to obtain the NDVI threshold in the alpine timberline and to determine the location of the alpine timberline. The Google Earth Engine platform provides common remote sensing image data sets that can be called and calculated on online, providing convenience for the long-term monitoring of alpine timberlines [42,43]. Given these problems and new method, we calculated the NDVI based on Landsat imagery from 1987-2019, applied the Otsu method to Google Earth Engine, and obtained climate raster data by spatial interpolation that was based on terrain and generated topographic data by spatial analysis in the ArcGIS 10.2 software (ESRI, Inc., Redlands, CA, USA) across Gongga Mountain. Our primary objectives were: (1) determining the NDVI threshold of alpine timberline via the Otsu algorithm, (2) analysing spatial and temporal dynamics of the alpine timberline and exploring the effects of topography on alpine timberline distribution, and (3) revealing the comprehensive relationships among alpine timberline, topographic and climatic factors and discussing the variations of alpine timberline in future.

Study Area
The highest peak of Gongga Mountain (29°22′-29°57′N, 101°45′-102°10′E) is 7556 m above sea level ( Figure 1). (c) Gongga Mountain location. Note: the "1", "2", "3", and "4" denote the selected hillsides (i.e., first, second, third, and fourth hillside). The first and the third slopes denote sunny slopes, the second slope denotes a shady slope, and the fourth hillside denotes a mixed slope. (c) Gongga Mountain location. Note: the "1", "2", "3", and "4" denote the selected hillsides (i.e., first, second, third, and fourth hillside). The first and the third slopes denote sunny slopes, the second slope denotes a shady slope, and the fourth hillside denotes a mixed slope. Gongga Mountain features a geomorphic high mountain and deep valley landscape with average annual temperature and precipitation of 4.2 • C and 1947 mm, respectively [44,45], and possesses a transitional seasonal climate between the subtropical humid monsoon and the cold Tibetan Plateau climates, which is dominated by the southeast Pacific monsoon [30,32]. Because of the complex and varied natural geographical conditions, ecosystem types and biodiversity are extremely rich on Gongga Mountain. The integrity of natural altitudinal zonation ranging from subtropical to alpine vegetation is unique, making Gongga Mountain become an ideal area for studying alpine timberline in China [46]. In addition, Gongga Mountain is a national nature reserve surrounded by scenic spots. Of these, the Hailuogou is a national forest park, as the national 5A level scenic spot, which is the only glacier forest park in China and is a very popular tourist destination.

Landsat Satellite Imagery
The long-term monitoring of vegetation and land cover changes is supported by Landsat satellite imagery, which provides uninterrupted global multispectral surface imagery in 16-d intervals from July 1972 to the present [47]. We calculated the NDVI from 1987-2019 using satellite images with a 30 m spatial resolution from Landsat 5, 7, and 8. We discarded images from 2002, 2010, and 2017 because of severe cloud cover and missing image.

Topographic Data
The Digital Elevation Model (DEM) data are collected from the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) Version 3 (V3) provided by NASA EARTH DATA (https://earthdata.nasa.gov) [48]. The spatial resolution of the ASTER GDEM image is 30 m. The DEM vertical accuracy and horizontal resolution are 12.1 and 72 m, respectively. The DEM data are processed with ARCGIS10.2 software to obtain topographic factors (i.e., slope and aspect).

Google Earth Engine
Google Earth Engine (https://earthengine.google.com)-a platform with massive computational capabilities that can solve many geospatial analysis [42]. In this cloud-based platform, multiple researchers can simultaneously access open remote sensing product data and perform geospatial analysis on planet scale [43]. Common remote sensing products are efficiently processed through massive parallelisation using JavaScript or Python. We calculated the NDVI from Landsat satellite imagery and extracted the alpine timberline by the code editor on Google Earth Engine.

Otsu Method
In image interpretation, we usually want a method to automatically differentiate relatively homogenous categories such as land and water or forest and grass [49]. For single band images, two types classification can be implemented by finding an optimal Otsu threshold [39][40][41]. The Otsu method is an automatic and efficient way to realise cluster-based automatic image segmentation without depending on other priori information; furthermore, Otsu method possesses the flexible adaptation for different data and allows to monitor different ground targets including herbaceous crops, forest and surface water [50,51]. In this study, we used the Otsu method to identify the alpine timberline. Optimal thresholds are calculated by maximising the inter-class variance (ICV) of the pixel value distribution. The average grey degree and the inter-class variance are defined as follows [39]: where µ and ICV express the average grey degree and inter-class variance, respectively; n 0 is the number of pixels less than the threshold; n 1 is the number of pixels more than the threshold; n is the total number of pixels in the image; µ 0 and µ 1 are the average greyscale of the target and non-target categories, respectively; and ω 0 and ω 1 are the proportions of target and non-target pixels in the total pixels, respectively. The optimal threshold is obtained by traversing each grey value and comparing the inter-class variance. We conducted Otsu analysis on vegetation using the NDVI, and the obtained threshold was used as a criterion for extracting the alpine timberline. In order to extract alpine timberline better, we calculated the NDVI for the early growing season with Landsat images. The NDVI is defined as: where NIR and Red are the spectral reflectance values calculated in the near-infrared and red bands, respectively. We screened out regions with NDVI greater than 0.1 in Landsat images for the Otsu method analysis to reduce the impact of non-vegetated areas [52,53] and ultimately acquired the optimal Otsu threshold.

Alpine Timberline Extraction
We smoothed the NDVI images by convolving with a Gaussian kernel (which is generated from a continuous Gaussian signal on Google Earth Engine) to obtain NDVI threshold pixels and to detect dynamic alpine timberline changes. The contour line was obtained by the edge detection of NDVI pixels. We used a zero-crossing edge detection method, which identifies the zero-crossings value on each frequency band of the image. Next, we selected the contour line with a property value equal to the threshold value; pixels on the contour line are the alpine timberline points. The alpine timberline points were also post-processed (discussed later in Section 2.4).

Slope and Aspect Classification
Gongga Mountain has a marked vertical span of more than 6000 m. Slope, a significant geomorphological parameter, is the angle between the land surface and the horizontal plane expressed in degrees. Different slopes may affect the growth of trees, and higher slopes have more unstable foundations. The slope on Gongga Mountain is categorised as flat, very gentle, gentle, moderate, moderately steep, steep, and very steep for slope ranges of 0-2 • , 2-5 • , 5-15 • , 15-25 • , 25-35 • , 35-55 • , and more than 55 • , respectively [54]. Table 1 presents the slope classification of Gongga Mountain. The compass direction on a topographic slope is definition of aspect, usually starts at geographic north. Aspect strongly influences the radiation intensity, sunshine hours, and temperature, water, and soil conditions of different slopes, affecting vegetation growth and the distribution characteristics of the alpine timberline. Flat ground has no slope or aspect.  [54]. Table 2 presents the aspect classification of Gongga Mountain.

Driving Force Analysis
The dynamics of alpine timberline are the result of the interaction of multiple factors. In this study, we chose four hillsides in the western part of Gongga Mountain and acquired the topographic and climatic data at the highest alpine timberline point from 1987-2019 to explore the relationship between alpine timberline, topographic and climatic factors. The west of Gongga Mountain was chosen for its relatively pristine environment conditions and suitable elevation (above 3000 m). We selected a total of 10 factors, including slope, aspect, longitude, latitude, average annual temperature, annual precipitation, growing season temperature, maximum growing season temperature, minimum growing season temperature, and growing season precipitation for each hillside. A combination of Pearson correlation analysis, principal component analysis, and multiple linear regression were performed.
(1) Pearson correlational analysis The Pearson correlation analysis is a classic mathematical method used to measure the degree of correlation between variables of two groups. The Pearson correlation coefficient quantifies this degree (ranging from −1 to 1) and is calculated as follows: where r is the Pearson correlation coefficient, x and y are the variables of the two groups, and n is the number of variables in each group. The variables of the two groups are more closely correlated with absolute values of r closer to 1 [55]. This analysis was performed using the "corrplot" package in the R software (R Core Development Team, R Foundation for Statistical Computing, Vienna, Austria).
(2) Interaction analysis of multiple factors The principal component analysis (PCA) is a common multivariate method that extracts important information from multidimensional data and visualises principal components in the R software [56]. The number of principal components is determined by the cumulative total variance contribution rate, which should be greater than 75%. Moreover, we used the two R packages "FactoMineR" and "factoextra" to analyse and visualise the data. Multiple linear regression (MLR) is a statistical technique that establishes a linear relationship between explanatory and response variables [57]. The results of the principal components were used as the explanatory variables of MLR in this study.

Processing
Data processing proceeded as follows ( Figure 2): first, we calculated the NDVI of Gongga Mountain from Landsat images in Google Earth Engine. Second, the optimal NDVI thresholds were calculated by the Otsu method, and the threshold pixels were extracted. Next, we limited the alpine timberline points to above 3200 m above sea level and verified the higher alpine timberline points to obtain the final alpine timberline points of Gongga Mountain based on historical Google images, land classification maps, and existing research results [26,32,58]. Then, the nearest alpine timberline points were obtained by seeking the alpine timberline points closest to the mountaintop point in the vicinity of each mountaintop point. After that, we analysed the spatiotemporal dynamic pattern of the nearest alpine timberline points and discussed the distribution of the nearest alpine forest line points on different terrains. Finally, we explored the relationships among the distribution of the highest alpine timberline points and topographic and climatic factors on four different hillsides from 1987-2019. absolute values of r closer to 1 [55]. This analysis was performed using the "corrplot" package in the R software (R Core Development Team, R Foundation for Statistical Computing, Vienna, Austria).
(2) Interaction analysis of multiple factors The principal component analysis (PCA) is a common multivariate method that extracts important information from multidimensional data and visualises principal components in the R software [56]. The number of principal components is determined by the cumulative total variance contribution rate, which should be greater than 75%. Moreover, we used the two R packages "FactoMineR" and "factoextra" to analyse and visualise the data. Multiple linear regression (MLR) is a statistical technique that establishes a linear relationship between explanatory and response variables [57]. The results of the principal components were used as the explanatory variables of MLR in this study.

Processing
Data processing proceeded as follows ( Figure 2): first, we calculated the NDVI of Gongga Mountain from Landsat images in Google Earth Engine. Second, the optimal NDVI thresholds were calculated by the Otsu method, and the threshold pixels were extracted. Next, we limited the alpine timberline points to above 3200 m above sea level and verified the higher alpine timberline points to obtain the final alpine timberline points of Gongga Mountain based on historical Google images, land classification maps, and existing research results [26,32,58]. Then, the nearest alpine timberline points were obtained by seeking the alpine timberline points closest to the mountaintop point in the vicinity of each mountaintop point. After that, we analysed the spatiotemporal dynamic pattern of the nearest alpine timberline points and discussed the distribution of the nearest alpine forest line points on different terrains. Finally, we explored the relationships among the distribution of the highest alpine timberline points and topographic and climatic factors on four different hillsides from 1987-2019.

Spatiotemporal Alpine Timberline Distribution
The height of the alpine timberline rises from 3203 to 4889 m around Gongga Mountain ( Figure  5a). There was spatial heterogeneity in dynamics of timberline. There had been a downward, upward, and unobvious trend in the northwest (Figure 5b

Spatiotemporal Alpine Timberline Distribution
The height of the alpine timberline rises from 3203 to 4889 m around Gongga Mountain ( Figure  5a). There was spatial heterogeneity in dynamics of timberline. There had been a downward, upward, and unobvious trend in the northwest (Figure 5b

Spatiotemporal Alpine Timberline Distribution
The height of the alpine timberline rises from 3203 to 4889 m around Gongga Mountain (Figure 5a). There was spatial heterogeneity in dynamics of timberline. There had been a downward, upward, and unobvious trend in the northwest (Figure 5b

Aspect Factor
The aspect distribution of the alpine timberline points is different on Gongga Mountain ( Figure  6). The northeast and west slopes have the highest (901) and lowest (177) Figure 8 shows the relationships between the elevation of the nearest alpine timberline points and the longitude and latitude. There was a significant relationship between the longitude and elevation of the alpine timberline (R 2 = 0.32, p < 0.01, Figure 8a). Conversely, the relation between latitude and elevation was not significant (binomial fitting: R 2 = 0.04, p < 0.01, Figure 8b). The elevation of the alpine timberline points generally drops gradually from west to east (Figure 8a), but the number of alpine timberline points increases from west to east. There are fewer alpine timberline points at high latitudes than at low latitudes (Figure 8b).

Slope Factor
We performed a statistical distribution of the nearest alpine timberline points on different slope classifications (Figure 7). The distribution of the nearest alpine timberline points is mostly on moderate    Figure 8 shows the relationships between the elevation of the nearest alpine timberline points and the longitude and latitude. There was a significant relationship between the longitude and elevation of the alpine timberline (R 2 = 0.32, p < 0.01, Figure 8a). Conversely, the relation between latitude and elevation was not significant (binomial fitting: R 2 = 0.04, p < 0.01, Figure 8b). The elevation of the alpine timberline points generally drops gradually from west to east (Figure 8a), but the number of alpine timberline points increases from west to east. There are fewer alpine timberline points at high latitudes than at low latitudes (Figure 8b).  Figure 8 shows the relationships between the elevation of the nearest alpine timberline points and the longitude and latitude. There was a significant relationship between the longitude and elevation of the alpine timberline (R 2 = 0.32, p < 0.01, Figure 8a). Conversely, the relation between latitude and elevation was not significant (binomial fitting: R 2 = 0.04, p < 0.01, Figure 8b). The elevation of the alpine timberline points generally drops gradually from west to east (Figure 8a), but the number of alpine timberline points increases from west to east. There are fewer alpine timberline points at high latitudes than at low latitudes (Figure 8b).

Driving Factors of the Spatiotemporal Distribution of Alpine Timberlines
On the first hillside, longitude had the highest correlation coefficient with elevation (r = 0.5, p < 0.01, Figure 9a).  Figure 1 shows the locations of hillsides 1-4; ** and * indicate significance coefficients of less than 0.01 and 0.05, respectively. According to Figure 10a,b, the highest R 2 of MLR was 0.122 based on two principal components (Dim2 and Dim3). Drivers with a greater contribution included longitude, latitude, Year_PRE, and

Driving Factors of the Spatiotemporal Distribution of Alpine Timberlines
On the first hillside, longitude had the highest correlation coefficient with elevation (r = 0.5, p < 0.01, Figure 9a).

Driving Factors of the Spatiotemporal Distribution of Alpine Timberlines
On the first hillside, longitude had the highest correlation coefficient with elevation (r = 0.5, p < 0.01, Figure 9a).  Figure 1 shows the locations of hillsides 1-4; ** and * indicate significance coefficients of less than 0.01 and 0.05, respectively. According to Figure 10a,b, the highest R 2 of MLR was 0.122 based on two principal components (Dim2 and Dim3). Drivers with a greater contribution included longitude, latitude, Year_PRE, and  Figure 1 shows the locations of hillsides 1-4; ** and * indicate significance coefficients of less than 0.01 and 0.05, respectively. According to Figure 10a,b, the highest R 2 of MLR was 0.122 based on two principal components (Dim2 and Dim3). Drivers with a greater contribution included longitude, latitude, Year_PRE, and Grow_PRE. Longitude of the second hillside possessed the highest r value with elevation (r = 0.49, p < 0.01), and the factor of the largest negative correlation was slope (r = -0.51, p < 0.01, Figure 9b). The interaction of Dim1, Dim3, and Dim4 contributed the highest R 2 value (R 2 = 0.152), which accounted for the Grow_TEM, Year_TEM, longitude, and Grow_minTEM (Figure 10c,d). For the third hillside, latitude (r = 0.74, p < 0.01), longitude (r = 0.55, p < 0.01), Grow_PRE (r = 0.52, p < 0.01), and Year_PRE (r = 0.49, p < 0.01) were the principal factors (Figure 9c). The highest R 2 value (R 2 = 0.49) was a consequence of the interaction of Dim1, Dim3, and Dim4; the principal contributed drivers were Year_PRE, Grow_PRE, latitude, and longitude (Figure 10e,f). The relationship between elevation and aspect factor was most significant on the fourth hillside (r = −0.44, p < 0.05, Figure 9d). The highest R 2 value occurred in Dim4, and the critical factors were Grow_maxTEM, slope, aspect, and Grow_minTEM (Figure 10g,h).

Temporal and Spatial Dynamics of Alpine Timberline on Gongga Mountain
The elevation of the nearest alpine timberline points to peaks on Gongga Mountain were extracted by remote sensing technology ranges from 3202-4889 m ( Figure 5) [48]. The extracted alpine timberline is slightly higher (or lower) than the previous research results on the alpine timberline of the southeastern margin of the Tibetan Plateau (3300-4600 m) [26,31,54,59], which could be attributed to the image resolution (30 m) and DEM vertical accuracy (±12 m). In the last 30 years, the NDVI threshold obtained by the Otsu method has changed significantly (Figure 4), indicating a significant increase in the vegetation coverage near the alpine timberline ecotone. The results are consistent with those of another study that found a significant increase in the population density of Abies fabri near the alpine timberline in the Yajiageng area of Gongga Mountain over the past 100 years [32]. Additionally, some studies showed that the population density of the Smith fir timberline in the southeastern Tibetan Plateau has significantly increased with climate warming since 200 years ago

Temporal and Spatial Dynamics of Alpine Timberline on Gongga Mountain
The elevation of the nearest alpine timberline points to peaks on Gongga Mountain were extracted by remote sensing technology ranges from 3202-4889 m ( Figure 5) [48]. The extracted alpine timberline is slightly higher (or lower) than the previous research results on the alpine timberline of the southeastern margin of the Tibetan Plateau (3300-4600 m) [26,31,54,59], which could be attributed to the image resolution (30 m) and DEM vertical accuracy (±12 m). In the last 30 years, the NDVI threshold obtained by the Otsu method has changed significantly (Figure 4), indicating a significant increase in the vegetation coverage near the alpine timberline ecotone. The results are consistent with those of another study that found a significant increase in the population density of Abies fabri near the alpine timberline in the Yajiageng area of Gongga Mountain over the past 100 years [32]. Additionally, some studies showed that the population density of the Smith fir timberline in the southeastern Tibetan Plateau has significantly increased with climate warming since 200 years ago [60]. This also happened in other areas, such as the Spanish Pyrenees and the central Canadian Rockies [61,62]. The harsh environment of the alpine timberline ecotone could affect the shape and size of vegetation but does not affect seed propagation and vegetation regeneration; this fact supports our results [31,60]. The alpine timberline generally shows a slow rising trend, but the rise is not obvious in some areas ( Figure 5), indicating the spatial heterogeneity of the alpine timberline. On a global scale, a meta-analysis of 166 datasets showed that 47% of treelines experienced no significant elevation changes [63]. For example, the elevation of alpine timberline did not significantly increase in the Yajiageng area because of growth restricted by grazing [32]. Furthermore, the advance of growing season and the increase of freezing events lead to alpine timberline stability [64,65].
The distribution of alpine timberline points differs on the different aspect. The elevation of the alpine timberline on sunny slopes is generally higher than that of shady slopes (Figure 6), because vegetation receives different solar radiation intensity and duration and heat availability depending on the aspect [66,67]. Additionally, the slope also affects the distribution of alpine timberline (Figure 7). On the one hand, the distribution of alpine timberline points on Gongga Mountain are relatively few in flat and very gentle places (Figure 7a)-because vegetation generally grows better on these slopes. Meanwhile, these flat and extremely gentle slopes are rare on Gongga Mountain (Table 1). On the other hand, fewer alpine timberline points are distributed on very steep slopes (>55 • ) (Figure 7a), mainly due to the large amount of gravel, scant soil moisture, and unfavorable seed preservation [68,69]. Thus, alpine timberline points are less likely to occur on excessively high or low slopes. The elevation of the alpine timberline decreased from west to east because of the topography of Gongga Mountain-which is higher in the west and lower in the east (Figures 1 and 8a). The latitude in this area had little influence on the distribution of alpine timberline points (Figure 8b), mainly due to the small latitudinal direction span and the complex terrain.

Drivers Influencing the Spatiotemporal Distribution of the Alpine Timberline
Our findings showed that the drivers of the alpine timberline present spatial heterogeneity on the west slope of Gongga Mountain (Figures 9 and 10). On the first hillside, the longitude has the most significant correlation with the alpine timberline elevation (p < 0.01, Figure 9a); longitude, latitude, Year_PRE, and Grow_PRE exceed the average contribution rate for driving factor interaction (Figure 10b), indicating that longitude is the primary driving factor in this region. However, the alpine timberline migration mechanism must be further observed and studied because of the small terrain span [26]. The contribution rate of Grow_TEM, Year_TEM, longitude, and Grow_minTEM exceeded the average value for driving factor interaction on the second hillside, and precipitation contributed less than average (Figure 10d); slope was the largest negative correlation factor here (p < 0.01, Figure 9b), with increased slopes making water preservation difficult and soil nutrient loss to rainfall more likely, thereby affecting alpine timberline migration [70]. Precipitation is the principal driving force on the third hillside (p < 0.01, Figures 9c and 10e,f). The following reasons explain why precipitation has a significant impact on the alpine timberline; first, precipitation defines the nutrient availability and growth restrictions of vegetation in alpine timberline ecotones [71][72][73]. Second, vegetation possess the limited capacity of precipitation use efficiency on the Tibetan Plateau, leading to wilted seedlings and limited tree growth at high altitudes [74][75][76]. The alpine timberline of the fourth hillside is directly impacted by aspect (p < 0.05, Figure 9d), suggesting that alpine timberlines are constrained by the intensity of solar radiation and sunshine hours, which is consistent with the research results of the timberline ecotone in the Changbai Mountains [77]. Moreover, previous study has shown that the vegetation distribution on Gonga Mountain is controlled by slope [78]. Interestingly, based on the contribution degree, we found that the factors exceeding the average contribution on the first and third hillsides are the same, whereas the factors exceeding the average contribution of the second hillside are completely different to those ( Figure 10). Because the first and the third slopes are sunny slopes and the second slope is a shady slope, the alpine timberline differs significantly among these slopes [32]. According to the analysis of these four hillsides, the drivers and migration of the alpine timberline have the hillside differences, and depend on climatic variation and the complex terrain, primarily different hydrothermal conditions, which affects the inconsistency of vegetation growth [79,80].
We conducted a literature survey addressing drivers of alpine timberline (Table 3). Regional differences of driving factors show dissimilarly dynamics of alpine timberline, but water and heat are the most basic factors affecting the alpine timberline. The accelerated retreat of glaciers caused by climate warming; this could lead to a more significant upward alpine timberline trend on Gongga Mountain [81,82]. However, species interactions may prevent the migration of the alpine timberline to higher elevations of the Tibetan Plateau [83]. Predicting and simulating alpine timberline dynamics is increasingly difficult because of the frequency of extreme weather, human activities and the lagging and diverse response of the alpine timberline to climate change [84][85][86]. Therefore, a combination of biotic and abiotic factors should be considered when predicting and simulating alpine timberline migration. Table 3. Alpine timberline driving factors in different regions.

Region
The Drivers Alps Maximum month temperature [13] Siberia Summer temperatures and growing season length [25] Winter temperatures and precipitation [27] Himalaya Growing season and annual bio-temperature [26] Topography and human disturbances [87] Temperature and moisture supply [88] Carpathians Growing season temperatures and late winter precipitation [89] Land use, grazing, logging, air pollution, and climate warming [90] Sierra Nevada Precipitation [91] Northwest Alaska Microtopography [92] Rocky Mountains Geomorphology and Geology [93] Global scale Winter warming [63] Low temperature, average temperature, and average soil temperature during the growing season [88]

Importance and Uncertainties
In this study, the Otsu method was used to extract the alpine timberline on the spatial scale that based on the Google Earth Engine and calculate the Otsu threshold of NDVI across the Gongga Mountain during 1987-2019, which captured alpine timberline variation on the pixel level and discussed the drivers of alpine timberline. The Otsu method is very adaptable to different time scales, and regions, and it can also confirm the optimal threshold of alpine timberline. For the extracting alpine timberline, we selected images derived from the early growing season, which can capture the growth characteristics of different vegetation types based on the differences of NDVI. Additionally, the ASTER GDEM V3, highly accurate DEM, was introduced to calculate the elevation of alpine timberline [48], which can improve the identification accuracy of alpine timberline. As we known, the climate is easily affected by the complicated topography, which is indispensable datasets for studying the drivers of alpine timberline; fortunately, the Anusplin 4.2 software is good tool which can improve the spatial interpolations accuracy in term of variations of topography. Moreover, Google Earth Engine provides better opportunities for undertaking alpine timberline observation studies because of its huge datasets and efficient computing power [42,43]. Ultimately, the drivers of different hillsides highlighted the different response processes and mechanisms of alpine timberline on spatial scale, which is helpful for future studies regarding the effect of climate change and topography on dynamics of alpine timberline.
Despite we obtained a certain achievement, the difference in the spectral bands and the instrument performances between Landsat 8 OLI and Landsat-5 TM/Landsat-7 ETM+ might cause inconsistency on the derived NDVI [94,95]. However, some studies have shown that NDVIs of Landsat 8 OLI and Landsat-7 ETM+ have good consistency, especially for vegetation types, they have better NDVI consistency than non-vegetation types [96,97]. Besides, it is significant challenges to extract alpine timberline on the pixel level with vegetation composites in thematic map because of the quality and the uneven availability of Landsat images at spatiotemporal scales [98,99]. Consequently, we checked multi-year alpine timberline maps by historical Google images, land classification maps, and existing research results, and minimized the uncertainties that derived from the data availability. Furthermore, vertical and horizontal errors in the DEM could led to the overall uncertainty for the alpine timberline elevation. In our study, the compounded dynamic mechanism study of timberline is scarce due to the limitation of data. Thus, we used previous research to explain our results, and the specific in situ observation should be considered in the future. If additional data becomes available, further studies will be conducted to evaluate the performance of the Landsat maps with alpine timberline elevation, density and future trend at the sub-pixel scale.

Conclusions
Alpine timberlines are warning signal and amplifier of global change. Remote sensing and geographic information system are the most essential tools for studying alpine timberline dynamics. In this research, we emphasized that the alpine timberline is an ecosystem ecotone rather than a vegetation type boundary. We used the Otsu method to extract the alpine timberline on Google Earth Engine, analysed the spatiotemporal pattern of the alpine timberline, and detected the controls of alpine timberline using Pearson's correlation coefficient, MLR, and PCA. Our results suggested that the vegetation coverage near the alpine timberline increases notable. The alpine timberline rose slowly with insignificant change in some areas and is controlled by aspect, slope, longitude, and latitude. The drives of alpine timberline had spatial heterogeneity on Gongga Mountain. Notable differences exist on sunny and shady slopes of alpine timberline. Consequently, more attention should be devoted to monitoring and simulating dynamics of alpine timberline on different hillsides. Additionally, the physical and chemical properties of soil and vegetation traits should be investigated in the field to better monitor the changed mechanisms of alpine timberline and to further deepen the understanding of the linkage of the alpine ecotone with climate change.