Detecting the Turning Points of Grassland Autumn Phenology on the Qinghai-Tibetan Plateau: Spatial Heterogeneity and Controls

: Autumn phenology, commonly represented by the end of season (EOS), is considered to be the most sensitive and crucial productivity indicator of alpine and cold grassland in the Qinghai-Tibetan Plateau. Previous studies typically assumed that the rates of EOS changes remain unchanged over long time periods. However, pixel-scale analysis indicates the existence of turning points and differing EOS change rates before and after these points. The spatial heterogeneity and controls of these turning points remain unclear. In this study, the EOS turning point changes are extracted and their controls are explored by integrating long time-series remote sensing images and piecewise regression methods. The results indicate that the EOS changed over time with a delay rate of 0.08 days/year during 1982–2015. The rates of change are not consistent over different time periods, which clearly highlights the existence of turning points. The results show that temperature contributed most strongly to the EOS changes, followed by precipitation and insolation. Furthermore, the turning points of climate, human activities (e.g., grazing, economic development), and their intersections are found to jointly control the EOS turning points. This study is the ﬁrst quantitative investigation into the spatial heterogeneity and controls of the EOS turning points on the Qinghai-Tibetan Plateau, and provides important insight into the growth mechanism of alpine and cold grassland.


Introduction
Vegetation phenology refers to periodically recurring growth patterns [1], and sheds a unique light on how ecosystems respond to climate change [2][3][4]. Shifts in phenology trends can affect the carbon budget, water flux, and energy balance from a regional to global scale [5]. Regional warming in alpine regions has led to several significant phenology changes, including advancement of the start of the growing season (SOS) in spring and a delay of the end of season (EOS) in autumn, as well as extensions of the growing season [6]. Phenology changes in turn provide strong feedback to climate systems, which can affect the regional carbon and water cycles [7]. The advancement of SOS and its riches the understanding of the EOS controls on alpine and cold grassland but also provides further details to reveal the EOS change mechanisms over different periods and their controls on the Qinghai-Tibetan Plateau.

Study Area
The Qinghai-Tibetan Plateau is situated in southwestern China and covers all of Tibet and the Qinghai provinces and is also a part of the Xinjiang, Sichuan, Gansu, and Yunnan provinces ( Figure 1). The Qinghai-Tibetan Plateau is considered the third pole in the world, has an average altitude of >4000 m, and is characterized by a plateau monsoon climate with low temperatures, low precipitation, and strong insolation. More than 54% of the Qinghai-Tibetan Plateau area has a total annual precipitation below 400 mm and temperatures below 0 °C [22]. This region is known as the Asia water tower and is home to the headstreams of the Yangtze, Yellow, Lantsang, and Indus rivers. The alpine, cold, and dry climatic conditions lead to unique vegetation types on the Qinghai-Tibetan Plateau. A climate gradient exists from warm-humid in the southeast to cold-dry in the northwest, along which the vegetation types transition from forestland, meadow, steppe, and desert. The grassland, which includes meadow, steppe, and desert steppe, and covers 51.05% of the Qinghai-Tibetan Plateau area, is the most important ecosystem and sensitive to climate change. An understanding of grassland dynamics under the climate and human disturbance conditions is crucial for regional ecological security. We divided the entire Qinghai-Tibetan Plateau into 12 subregions ( Figure 1, Table 1) based on the bio-geographical division proposed by Zheng et al. [30]. The grassland distribution was extracted according to a China vegetation map (scale = 1:100,000) [31], eliminating subregions X, XI, and XII, for which the main vegetation types are desert, forestland, and forestland, respectively. Only the remaining nine subregions (I-IX) are analyzed in this study, covering meadow, steppe, and desert grassland (Table 1). Of these nine subregions, we focused in detail on subregion I, which has the highest annual accumulated temperature above 0 °C (AGDD0) and medium moisture index (MI). Subregions We divided the entire Qinghai-Tibetan Plateau into 12 subregions ( Figure 1, Table 1) based on the bio-geographical division proposed by Zheng et al. [30]. The grassland distribution was extracted according to a China vegetation map (scale = 1:100,000) [31], eliminating subregions X, XI, and XII, for which the main vegetation types are desert, forestland, and forestland, respectively. Only the remaining nine subregions (I-IX) are analyzed in this study, covering meadow, steppe, and desert grassland (Table 1). Of these nine subregions, we focused in detail on subregion I, which has the highest annual accumulated temperature above 0 • C (AGDD 0 ) and medium moisture index (MI). Subregions II and III had relatively high MI values that decreased from southeast to northwest. Each subregion exhibited unique climatic conditions and economic development levels, as well as different vegetation responses to climate and human activities.

Data Source
The GIMMS NDVI3g dataset provided by NASA was used to estimate the EOS on the Qinghai-Tibetan Plateau. The dataset was available from 1982 to 2015 with an 8-km spatial resolution and 15-day temporal resolution [32]. Some previous processes (e.g., calibration, noise removal) were performed for this version to better detect the vegetation dynamics [32]. This dataset has been widely used to detect long-term vegetation dynamics [33][34][35]. Due to the normalized difference vegetation index (NDVI) data might be misrepresented by snow [36]; we used the average temperature of a sequence of five days less than 0 • C to screen out the pixels that might be covered by snow. Temperature, precipitation, and insolation data from 1982-2015 were extracted from the China meteorological forcing dataset  downloaded from the Big Earth Data Platform for Three Poles with a spatial resolution of 0.1 • and temporal resolution of 3 h (http://poles.tpdc.ac.cn/, Accessed on 15 August 2021) [37].
Human activities, including grazing density and economic development, were quantified using economic statistic data. For example, the grazing density were represented by the number of large animals (one large animal equal to five sheep unit) and sheep, and uniformly converted into sheep units. The economic development levels were quantified as the production of primary, secondary, and tertiary sectors. These data come from the statistical yearbooks of Qinghai and Tibet from 1982 to 2015.

Retrieval of EOS
Numerous methods have been used to fit the NDVI changes from seasonal vegetation cycles. After comparing the fitting results of HANTS [38], Polyfit [39], and Double logistic [40] in the nine subregions, we found that the RMES of HANTS (1.26 ± 0.24 × 10 −5 ) and Polyfit (1.28 ± 0.24 × 10 −5 ) were similar and smaller than the Double logistic results (1.93 ± 0.43 × 10 −5 ) ( Figure A1). HANTS and Polyfit, were therefore selected as the two most simple and effective methods to fit the NDVI change curves. Dynamic thresholds were adopted to determine the EOS. Further details of these two fitting methods are described below.
The HANTS method involves the harmonic analysis of a time series, is adapted from the fast Fourier transform, and eliminates cloud noise using the least square method [38,40]. The HANTS methods can quickly smooth the data, remove outliers, and fill gaps of missing data. The following Equation (1) was used to fit the NDVI seasonal fluctuation curve: where t is the Julian date, a 0 is the average of all NDVI observations, and ϕ i and a i are the phase and amplitude of the curve, respectively. The Polyfit method adopts a polynomial function to fit the NDVI records [39]. The following sixth order Equation (2) is used to describe the NDVI curve: NDV I(t) = a 0 + a 1 t + a 2 t 2 + . . . + a 6 t 6 (2) where a 0 -a 6 are regression coefficients determined using the Levenberg-Marquardt method. The EOS was determined by the day when the smoothed curve of the 34-year mean passed a designated threshold. We first fitted the NDVI changes with HANTS and Polyfit methods and then calculated the NDVI ratio (described in Equation (3)) for 365 days with multi-year mean NDVI values, next detected the time t with the minimum NDVI ratio and used the corresponding NDVI(t + 1) at time (t + 1) as the NDVI threshold for the EOS. Finally, we obtained the EOS for 34 years using the threshold:

Quantification of the EOS Trends, Turning Points, and Controls
After extracting the EOS at the pixel scale, we first quantified the tendency of the EOS changes using greenness changes methods, and then detected the turning points using the piecewise regression method. The mean EOS values and EOS trends within the subregions were calculated as the EOS and EOS trends at the subregion level. The turning points at a subregion and province level were calculated by the majority values.
The EOS trends were calculated using the greenness rate of change [41]. The EOS was considered delayed if the slope was a positive value; otherwise, the EOS advanced.
where i is the order of the year, n is the number of years, NDVI i is the NDVI in the i th year, and the slope is the vegetation change rate. Alternatively, we can use the unary linear regression, in which the P values and confidence levels can be calculated. Turning points were identified by piecewise regression [42] analysis, as defined in Equation (5), which can be used to detect sudden and sharp changes in directionality. This method has been widely applied for analyzing vegetation dynamics [19,43,44].
where t is the order of the year, α is the estimated turning point of the vegetation change trend determined using the least square error method, β 1 and (β 1 + β 2 ) are the change rates before and after the turning points, respectively, and ε is the residual error. We performed t-tests to check the significance of the piecewise regressions. Redundancy analysis (RDA) is a powerful analysis technique that could be applied in separating the contributions of climate, human activities, and their intersections to the EOS changes. RDA is a method to extract and summarize the variation in a set of response variables that can be explained by a set of explanatory variables [45]. In this study, RDA was performed with the vegan package in R language [46]. In RDA, climatic variables or human activity variables were chosen as predictors to maximize the extent of their correlation with the EOS changes as the response variable. RDA had been widely used in ecology-related studies [47,48]. The turning points of human activities were also calculated with Equation (5). The relationships between the turning points of the EOS and climatic variables were quantified using partial regression analysis or the correlation coefficient.

EOS Spatial Distribution and Variation Characteristics
The obtained EOS presented high spatial heterogeneity across the grassland of the Qinghai-Tibetan Plateau during 1982-2015 ( Figure 2). The EOS results extracted using the HANTS and Polyfit methods were not consistent, but their spatial distribution trends were similar (Figure 2a,b). The mean multi-year EOS began on the 291th day of the year (end of September) and spanned nearly one month from the southeast to the northwest (Figures 2 and 3a). The EOS started early (around the 277th day) in subregion IX, which has the highest elevation and lowest AGDD 0 values, and started late (around the 300th day) in subregions II and III, which are characterized by relatively warm-humid conditions. In the central Qinghai-Tibetan Plateau (subregion V), the EOS occurred on the 295th day. In subregion I, the EOS was early in the west and late in the east with a mean EOS on the 292th day. The spatial heterogeneity variations were significantly controlled by the MI (EOS = 16.55 × MI + 287.28, R 2 adj = 0.20 ** ), with an early EOS in the drought subregions (IV, VII, VIII, and IX) and late EOS in the relatively humid subregions (II, III, and VI). The EOS spatial heterogeneity was essentially insensitive to AGDD 0 .
The mean EOS on the Qinghai-Tibetan Plateau exhibited a slow delayed trend with an average rate of 0.08 days/year. The EOS results extracted using the HANTS and Polyfit methods presented similar patterns (Figure 2d,e). Using these two fitting methods, 60.2% of the study area presented delay trends (27.8% area is significant), while 39.8% of the study area presented advance trends (13.4% area is significant). The EOS trends differed between the nine subregions during 1982-2015 (Figure 3b), showing a delay in the northwest and an advance in the southeast. Subregions I and IX showed significantly delayed trends with more than 0.20 days/year. The EOS of subregion II, with a main land use type of wetland, was also delayed by a rate of 0.12 days/year. The EOS in subregion VIII, which is characterized by alpine, cold, and dry climatic conditions, presented a negative trend with the fastest variation rate (−0.12 days/year) compared with the other subregions. The EOS of subregion IV showed an advanced trend in the north but delayed trend in the south, with a mean EOS trend of 0.02 days/year. The EOS in subregions III and V showed slight advanced trends of -0.02 and −0.01 days/year, respectively. Subregions VI and VII both presented a slightly delayed trend with an average of approximately 0.04 days/year.

Detection of EOS Turning Points in the Subregions
The EOS changed over time and presented delayed trends during 1982-2015, but the rates of change were not fixed in each subregion over different periods, and notable turning points were observed (48.2% is significant) (Figure 4c). For example, in subregion I, the turning point occurs in the year 1994, for which the EOS was delayed before 1994 and slightly advanced after 1994 (Figure 4f). Similarly, subregion II showed a delayed EOS before 2002 and then a slightly advanced EOS after 2002. In the remaining subregions (III, IV, VI, and IX), the change trends were similar and the turning point year was 1994, where the EOS was delayed prior to 1994, suddenly advanced in 1995, and then maintained the previous change trend until 2015. The turning point trends in subregions V, VII, and VIII occurred in 1994, 1994, and 1999 respectively, but were not significant. These results demonstrate that the EOS changes clearly exhibit turning points and a wide range of EOS change trends with significant spatial heterogeneity on the Qinghai-Tibetan Plateau. The pattern of EOS turning points extracted by HANTS and Polyfit (Figure 4a,b,d,e)       previous change trend until 2015. The turning point trends in subregions V, VII, and VIII occurred in 1994, 1994, and 1999 respectively, but were not significant. These results demonstrate that the EOS changes clearly exhibit turning points and a wide range of EOS change trends with significant spatial heterogeneity on the Qinghai-Tibetan Plateau. The pattern of EOS turning points extracted by HANTS and Polyfit (Figure 4a,b,d, and e) have a small difference in subregion I and VI.

EOS Variations Controlled by Climatic Variables before and after Turning Points
The EOS changes exhibited close relationships with the climatic variables, but the dominant climatic variable differed in each subregion before and after its associated turning point ( Figure 5). Temperature was the dominant control over the EOS changes in most subregions (I, II, IV, VI, VII, VIII, and IX) before and after the turning point year. In contrast, subregion III showed that the EOS was mainly controlled by the precipitation. Central subregion VII showed that the EOS was jointly controlled by the effects of temperature and precipitation. The area where the EOS changes was controlled by temperature covered the largest proportion, followed by precipitation and insolation (Figure 5d). The results indicate that the proportions controlled by each climate variable changed before and after the turning point years. For example, the EOS in subregion V was controlled by precipitation before the turning point, which then switched to temperature ( Figure A2). The EOS of only approximately 40% of the area in subregion VI was significantly controlled by temperature prior to the turning point, which thereafter increased to 70%. Contribution of climates to EOS variation are similar with the fitting results of HANTs and Polyfit ( Figure A3). after the turning point years. For example, the EOS in subregion V was controlled by precipitation before the turning point, which then switched to temperature ( Figure A2). The EOS of only approximately 40% of the area in subregion VI was significantly controlled by temperature prior to the turning point, which thereafter increased to 70%. Contribution of climates to EOS variation are similar with the fitting results of HANTs and Polyfit (Figure A3).

Controls on the EOS turning Points
The changes of the annual EOS turning points are consistent with the turning points of the climate variables in most subregions ( Table 2). In subregions I and II, the year of the EOS turning point coincides with the year of the insolation and precipitation turning points, respectively. Furthermore, the years of the EOS and temperature turning points are consistent in subregions III-IX. The major determining climatic variable for the EOS turning points is precipitation, followed by temperature and insolation. The relationship with the EOS turning point and insolation is generally weak (R 2 < 0.05).

Controls on the EOS Turning Points
The changes of the annual EOS turning points are consistent with the turning points of the climate variables in most subregions ( Table 2). In subregions I and II, the year of the EOS turning point coincides with the year of the insolation and precipitation turning points, respectively. Furthermore, the years of the EOS and temperature turning points are consistent in subregions III-IX. The major determining climatic variable for the EOS turning points is precipitation, followed by temperature and insolation. The relationship with the EOS turning point and insolation is generally weak (R 2 < 0.05). The relationship between the EOS and human activities was studied at the province level owing to limited statistical data in certain counties and subregions. The economic data show a consistent turning point with the EOS. Before the turning point year (~1996 for Qinghai and~1995 for Tibet), Qinghai maintained a large amount of sheep, which reflected high grazing activity, and the economic development was slow with low production in the primary, secondary, and tertiary sectors. However, after the turning point year, the grazing intensity decreased and reached a stable change rate, whereas the economy developed rapidly, especially in the secondary sector. For Tibet, the grazing intensity was small before the turning point year but showed a rapid growth rate after the turning point. Similar to Qinghai, Tibet experienced fast economic growth after the turning point, especially in the tertiary sector.
At the province level, the annual EOS was closely related to climate, human activities, and their intersections (Table 3). For Qinghai, a combination of the turning points of climate and human activities can explain 78.86% of the EOS turning points changes, with climate independently accounting for 40.22% and human activities accounting for 10.45%. The intersections of climate and human activities can explain 28.19% of the EOS variation in Qinghai. In Tibet, the EOS change due to climate (66.17%) was larger than that in Qinghai and the effect of human activities (6.8%) was weaker. The climate and human activities intersections in Tibet (9.98%) were also smaller than in Qinghai.

Controls on the EOS and EOS Turning Points
This study is the first to demonstrate pixel-scale spatial heterogeneity of the EOS turning points and explain the turning point controls. The results indicate that the joint effects of climate variables and human activities are the main controls of the EOS turning points. The response of the EOS to environmental changes is complex. Some previous studies indicated that temperature plays a crucial role in EOS regulation [49] however, we show that the temperature control over the EOS is regulated by precipitation and insolation in the meadow and grassland ecotones. The cause of the turning points in most subregions is the abrupt change of temperature and precipitation. The results also reveal that insolation contributed considerably to the EOS changes, which is consistent with some previous reports that the EOS and its relation with GPP is mainly limited by insolation [50]. Other studies have reported that meadow shrinkage, decreased land cover, land albedo changes, and permafrost and seasonal frozen soil dynamics intersect with climate change, which alters the EOS trends [51].
Grazing is the most important human activity that affects grassland dynamics on the Qinghai-Tibetan Plateau [52]. The spatial heterogeneity of community increases, community function alteration, and biodiversity loss are considered to be some of the key disturbances that result in grassland degradation [53,54]. The pika population could also increase the effects of animal distribution on vegetation [55]. Overgrazing reduces the vegetation biomass and height, and restricts the regrowth ability of grassland. Our analysis shows that grazing activities in Qinghai notably decreased around 1998, coinciding with the implementation of national conservation policies (e.g., ecological compensation, restoration of degraded grassland). Grazing in Tibet was not active before 1995 and then rapidly increased, however grazing decreased after 2005 due to the late implementation of ecological conservation projects. The primary industry (mostly agriculture and animal husbandry) increased by nearly a factor of five in 1996-2015 compared with that in 1982-1994, which is also consistent with the EOS change turning points. The tertiary industry in Qinghai and Tibet quickly increased after the turning points, which indirectly reflects the intensification of human activities on the Qinghai-Tibetan Plateau.

Ecological Significance of the EOS and Its Turning Points
Phenological changes have great effects on the structure and function of ecosystems. At the community level, various species have different phenological responses to climate change, whereas the EOS can lead to a change in the competition for light and water conditions [17,56]. Moreover, plant species changes in the community introduced by the EOS can lead to phenological mismatches; for example, the period of high consumer demand for a resource does not match with the period of resource abundance [57]. At the ecosystem level, phenological grassland changes can modify certain land surface parameters (e.g., albedo, sensible heat flux, evaporation, boundary layer conductivity), which affects the regional carbon and water cycles [58]. For example, a later EOS may promote GPP and cause plants to close their stomates and increase water use efficiency if a soil deficit exists [59]. Moreover, the delayed EOS may also increase transpiration and partly offset the GPP, therefore leading to closer relationships between the net ecosystem productivities and EOS changes [60].
The existence of turning points indicates that the EOS trend over long-time periods does not remain unchanged, and the rates of EOS changes differ before and after these points. This observation has several advantages in ecosystem-related studies. First, climatic controls on the EOS in the Qinghai-Tibetan Plateau intersect with each other and follow non-linear relationships with the EOS. An analysis of the EOS before and after the turning points therefore helpful to evaluate the climatic driving mechanisms of the EOS. Second, the detection of spatial heterogeneity of the turning points is helpful for evaluating the large-scale implementation effects of ecological conservation projects. Third, an analysis of the turning points of the EOS relationships with ecosystem functions and services provide important guidelines for fine ecology planning and the development of protection policies.

Uncertainties, Challenges, and Future Directions
The uncertainties in this study arise from three aspects. First, although the EOS trends are consistent with the findings of MODIS NDVI and SPOT NDVI, some design shortcomings in the AVHRR sensor may potentially introduce noise into the GIMMS 3g NDVI dataset. Second, the human activities are difficult to quantify for lack of grazing data (intensity and boundary) and statistic data on the county levels for a long time. Third, there is a limited number of phenological stations on the Qinghai-Tibetan Plateau, and most are distributed in the east, which thus does not represent the EOS changes of the entire plateau. The results of the EOS extraction are not fully calibrated by observations owing to limited data availability.
We recommended the following perspectives for future studies. First, extreme climate events (e.g., cold, frost, drought) may have a more direct effect on vegetation phenology than gradual changes in mean climatic conditions [27,28]. Non-structural carbohydrate storage in plants is helpful to avoid damage caused by extreme events [61]. However, extreme climate conditions with variable frequencies and intensities in different seasons on the Qinghai-Tibetan Plateau require rigorous quantification. Second, although many studies have quantified the effects of climate variables in different seasons, spring phenology, growth season length, and human disturbances on the EOS changes, the joint contribution of these variables is low and the control mechanisms of the EOS and its turning points remain poorly understood. The strengthening and development of phenological observations stations are therefore necessary to explain the mechanism of phenology changes in the Qinghai-Tibetan Plateau. Third, ecosystem models are essential tools for simulating the carbon cycle in both historic and future climate scenarios however, their accuracies remain limited by the understanding of the EOS [62]. More reasonable algorithms and reliable observations are required to calibrate the ecosystem models, which will ultimately provide a new research direction but presently faces serious challenges.

Conclusions
This study applied multiple statistical methods and long-time series remote sensing data to determine the spatial heterogeneity and controls of autumn phenology on the Qinghai-Tibetan Plateau. The results are summarized as follows. (1) EOS turning points exhibit notable spatial heterogeneities. (2) The climatic controls of the EOS before and after the turning points varied in different subregions on the Qinghai-Tibetan Plateau. (3) Changes in the turning points are controlled by the joint effects of climate and human activities (grazing and economic development). This study is the first to demonstrate the spatial heterogeneity of turning points at a pixel scale and discuss their controls on the Qinghai-Tibetan Plateau, which is useful for exploring the mechanism of EOS changes and developing regional ecosystem conservation measures.   Figure A1. The NDVI fitting results of three equations, including HANTS (red lines and numbers), Polyfit (green lines and numbers), and Double logistic (red lines and numbers). (a-i) represent subregions I -IX, respectively. The upper-left panels are corresponding NDVIratio changes and EOS with HANTS and Polyfit fitting methods. The upper-right numbers are RMESs for three fitting methods. Figure A2. Relative proportions of the three climate variables that contributed to the EOS in each subregion. The legend coloring is same as in Figure 5. The white portions indicate that the proportion is not significant (p > 0.05).