Dynamic Threshold of Carbon Phenology in Two Cold Temperate Grasslands in China

Plant phenology, especially the timing of the start and the end of the vegetation growing season (SOS and EOS), plays a major role in grassland ecosystem carbon cycles. As the second-largest grassland country in the world, China’s grasslands are mainly distributed in the northern cold temperate climate zone. The accuracies and relations of plant phenology estimations from multialgorithms and data resources are poorly understood. Here, we investigated vegetation phenology in two typical cold temperate grasslands, Haibei (HB) and Inner Mongolia (NM) grasslands, in China from 2001 to 2017. Compared to ground vegetation phenology observations, we analyzed the performance of the moderate resolution imaging spectroradiometer MODIS phenology products (MCD12Q2) and two remote sensing-based vegetation phenology algorithms from the normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) time series (five satellite-based phenology algorithms). The optimal algorithm was used to compare with eddy covariance (EC)-based carbon phenology, and to calculate the thresholds of carbon phenology periods (SOSt and EOSt) in each site. Results showed that satellite-based phenology estimations (all five algorithms in this study) were strongly coupled with the temporal variation of the observed phenological period but significantly overestimated the SOS, predicting it to be over 21 days later than the field data. The carbon phenology thresholds of HB grassland (HB_SOSt and HB_EOSt) had a significant upward trend, with the multiyear average values being 0.14 and 0.29, respectively. In contrast, the thresholds of NM grasslands (NM_SOSt and NM_EOSt) also showed a certain upward trend, but it was not significant (p > 0.05), with the multiyear average values being 0.17 and 0.2, respectively. Our study suggested the thresholds of carbon phenology periods (SOSt and EOSt, %) could be simply and effectively estimated based on their significant relationship with the EC-based maximum of gross primary productivity observations (GPPmax) at a specific site and time. Therefore, this study suggested the thresholds of carbon phenology were not fixed even in a specific ecosystem, which also provided simple bridges between satellite-based vegetation phenology and EC-based carbon phenology in similar grasslands.


Introduction
The phenological phase of plant growth, especially the timing of the start and the end of the vegetation growing season (SOS and EOS), plays a key role in terrestrial ecosystem carbon and nutrient cycles [1][2][3][4][5]. Previous studies across broad disciplines have observed phenology shifting both in individual plants and land surfaces [6]. Many studies have documented close relationships between plant phenology phases and climate changes [6,7]. Specifically, in high-latitude and colder regions, the temperature acts as a greater trigger on plant phenology due to changes as a result of the necessary heat requirements of vegetation [8,9]. It is well known that warming advances leaf unfolding and delays leaf coloring [10], but this trend is generally spatially-heterogeneous [11][12][13]. Besides, the phenology response to climate change also shows unfixed trends [14,15], which have decelerated or even reversed in recent years [10,14,16]. In addition, plant phenology is not only affected by temperature but is impacted by many confounding factors (i.e., water stress, sunshine, nutrients, and snow cover, etc.) [7,[17][18][19][20]. Consequently, the quantification of such complex impacts from multiple drivers on plant phenology changes remains challenging [14]. Therefore, diverse observation tools and phenology detection methods from multiaccess data resources are needed to improve the accuracy of the plant phenology studies.
Remote sensing-based vegetation phenology detection from temporal series of multiple vegetation indices has supported major advances in plant phenology research [14,21]. The normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) are the two most used vegetation indices for phenological studies [12,13,15,22]. Generally, the approaches evaluating the land surface phenology are of two types, threshold-based and the relative change rate algorithm [23]. Threshold-based vegetation phenology algorithm monitors the phenological transition periods (SOS and EOS) through different threshold settings (a fixed or a relative value derived from the maximum and minimum) of the vegetation index at different spatial-temporal scales. For example, a fixed NDVI threshold of 0.45 and ecosystem type-dependent values (from 0.11 to 0. 35) were respectively used to detect the landscape phenology in eastern US and China [24,25]. In contrast, half (or another quintile) way between the minimum and maximum NDVI was defined to track the SOS and EOS in different ecosystems across the globe [10,[26][27][28]. The relative change rate algorithm used the temporal trends of the vegetation growth curves to identify the SOS and EOS [23,29]. In recent years, relative change rate algorithms were more widely applied, using logistic [26,[29][30][31], piecewise regression [10], and polynomial algorithms [8,12,25,32]. However, the performance and accuracy of multialgorithms is not well known, especially not for scientific comparisons among similar ecosystems.
Ground observed phenology from specific species is a direct and important way to validate remote sensing-based vegetation phenology detection [33,34]. Eddy covariance (EC) technique, the most efficient auto-micrometeorological method to observe instantaneous and seasonal net CO2 exchange (NEE) between the biosphere and the atmosphere [35,36], is increasingly used for ecosystem model calibration and validation [37][38][39][40]. ECbased carbon observations provide another data source for tracking canopy phenology whose growing season is considered as the number of days when the ecosystem is a net carbon sink (NEE < 0) (namely "carbon uptake phenology") [2,41,42]. However, it is difficult to accurately quantify vegetation SOS and EOS using carbon phenology (NEE < 0, i.e., carbon uptake period) because of the mismatch between vegetation growth curves (i.e., leaf development, flowering, and fruiting) and seasonal NEE. Thus, many studies used threshold-based or inflexion-based approaches to track the EOS and SOS in seasonal series of gross primary productivity (GPP) derived from EC-based observations [26,43]. However, the threshold of GPP-based phenology estimation was variable across diverse ecosystems and at different spatiotemporal scales, which directly influenced the accuracy of phenology estimation.
The grassland ecosystem in China is approximately 400 million hectares, accounting for more than 40% of the land area, making China the second-largest grassland country across the globe [44,45]. China's grassland ecosystems are mainly distributed in the northern temperate continental arid climate zone and the western Qinghai-Tibet Plateau [44,46]. Therefore, in this study, we investigated vegetation phenology in two cold temperate grasslands in China. One was a typical grassland in the northern temperate continental arid climate zone (Inner Mongolia), and the other was in the Qinghai-Tibet Plateau (Haibei). Using two remote sensing-based algorithms and ground-based observations, the main objectives of our study were: (1) to evaluate the accuracy of remote sensing-based phenology estimation in temporal patterns and absolute amounts, and (2) to identify the threshold of GPP-based phenology estimation in these two cold temperate grasslands in China.

Site Description
The experimental sites were two typical grassland ecosystems in China: Haibei (HB) and Inner Mongolia (NM), which were near the local grassland ecosystem monitoring stations. The HB Grassland Station is in Xihai county (36°57'N, 100°51'°E, elevation 3140 m), Qinghai Autonomous Region (Figure 1). The local area is categorized as plateau monsoon climate with characteristics of strong radiation, low air temperature, and short cool summers. Climatic records from the past five decades show that the average annual air temperature is −0.1 C, with a coldest monthly mean of −14.1 C in January and highs of 11.5 C in July. Average annual precipitation is 391.9 mm, over 80% of which is concentrated between May and September. The dominant species of this alpine meadow grassland are Stipa krylovii, Kobresia humilis, Cleistogenes squarrosa, and Artemisia frigida, and the canopy heights are 20-30 cm. The soil is chestnut soil. The NM grassland station in Xilinhot, a representative site of the temperate steppe area, is in the Xilin Gol League in the middle of the Inner Mongolia Autonomous Region (43°57'N, 116°07'°E, elevation 1003 m) ( Figure 1). The climate is semiarid temperate continental climate. The coldest monthly mean of temperature is −18.1 C in January, and the warmest is 22.5 C in July. The average annual air temperature is higher than Haibei (2.5 C), whose ≥0 ℃ biological accumulated temperature can reach 2350～3400 days per year, but average annual precipitation is lower (273.7 mm). The dominant species of this temperate steppe are Stipa krylovii and Kobresia pygmaea, and the associated species are Leymus chinensis, Cleistogenes squarrosa, Agropyron cristatum, and Artemisia frigida. The soils are sandy loams.  [35]. EC flux observations in HB were from 2003 to 2011, and that in NM were from 2004 to 2011, respectively. To avoid potential influences from instrument malfunction, rainfall, dew or cobwebs, human disturbance, and near-static atmospheric conditions, raw measurements were processed and filled according to ChinaFLUX data processing [47,48], which included despiking, coordinate rotation, air density corrections, outlier rejection, and friction velocity threshold (u*) corrections [49]. Based on flux measurements and the exponential relationship between ecosystem respiration (Re) and soil temperature, daytime NEE was partitioned into EC-based GPP as carbon assimilation and Re as carbon emission during daytime (NEE = GPP + Re) [50,51].

EC and Meteorological Measurements
Standard meteorological and soil parameters were measured using an array of sensors, including photosynthetically active radiation (PAR), air temperature at the height of 2 m (Ta), precipitation (PPT), air relative humidity (RHa), vapor pressure deficit (VPD), soil temperature at depths of 5 cm (Ts), and soil water content (SWC). All measurements were averaged and recorded in a datalogger once per 30 min. The data were retrieved by a laptop computer every three weeks.

Phenology Measurements
The phenological monitoring of the dominant vegetation was conducted at the two sites, and the average value of the phenological period (SOS and EOS) of the dominant species in each site was taken as the phenological period of the grassland ecosystem. A 100 m × 100 m vegetation phenological observation field was selected and fenced according to the local grassland types in HB and NM, respectively. The observation field was usually divided into 4 observation areas (50 m × 50 m) as 4 repetitions. Four samples (1 m × 1 m) were randomly selected from each observation area. We selected 10 plants of each dominant species that were growing well and had a complete life history for 3 consecutive years as phenological observation objects. We considered the day of the year (DOY) that over 50% of each species (>5 plants) reached the SOS to the EOS as the phenology periods for the species. The phenological observation was conducted every 2 days for individual plants from SOS to EOS. The phenological observation periods of HB (HB_SOS and HB_EOS) and NM (NM_SOS and NM_EOS) were from 2001 to 2007 and from 2011 to 2017, respectively.

Remote Data Products
The moderate resolution imaging spectroradiometer (MODIS) Collection 6.0 land products are available from the Land Processes Distributed Active Archive Center (http://daac.ornl.gov/MODIS). In this study, using two site positions as the center pixel, surface vegetation index products (MOD13Q1), including NDVI and EVI, with a temporal resolution of 16 days and a spatial resolution of 250 m in the last 17 years from January 2001 to December 2017, were extracted. In addition, the surface phenology simulation products (MCD12Q2) from 2001 to 2017 were extracted.

Threshold Algorithm
In this study, based on the temporal dynamics of NDVI and EVI data in two research sites from 2001 to 2017, we used three methods, including asymmetric Gaussian function (GS), double logistic curve (LG), and Savitzky-Golay (SG) filtering, to fit and reconstruct the smoothing continuous time series of the two stations in TIMESAT software package [52]. According to the smoothed NDVI and EVI time series, the threshold method of TIMESAT software was used to extract the start time (SOS) and end time (EOS) of the growing period in the two stations, respectively. The threshold method is commonly used to simulate vegetation growth and can be expressed as Equation (1): where g(t; a1,a2…a5) represents the smoothed NDVI or EVI time series, a1 is the position (day of the year) of the minimum or maximum value of the smoothed temporal NDVI or EVI series, a2 and a3 are the width and flatness of the right half of the series, and a4 and a5 are the width and flatness of the left half of the series. Thus, the phenology estimation of the threshold method can be calculated by the average value of phenology calculated by three filtering method extractions (Figure 2a,c,e,g). The phenological results (SOS and EOS) calculated in HB and NM based on different vegetation indices (NDVI and EVI) are marked as HB_SOS_tn, HB_EOS_tn, HB_SOS_te, HB_EOS_te, and NM_SOS_tn, NM_EOS_tn, NM_SOS_te, NM_EOS_te, respectively (Table 1).

SOSt and EOSt
The vegetation carbon phenological period thresholds (%) of SOS and EOS.

Relative Change Rate Algorithm
In this study, the piecewise logistic function method was used to fit the time series of vegetation indices at two sites [29]. Specifically, the time change of vegetation index in a single growth or senescence cycle can be simulated using Equation (2): where y(t) is the NDVI or EVI value at the time of t, a and b are fitting parameters, c + d is the maximum value of NDVI and EVI in the year, and d is the initial value of NDVI and EVI (that is, the minimum value in the year). Among them, the maximum and minimum of the relative change rate of vegetation index are the beginning and ending points of the phenological period (Figure 2b,d,f,h). The phenological results (SOS and EOS) calculated by HB and NM based on different vegetation indices (NDVI and EVI) are marked as HB_SOS_rn, HB_EOS_rn, HB_SOS_re, HB_EOS_re and NM_SOS_rn, NM_EOS_rn, NM_SOS_re, NM_EOS_re, respectively (Table 1).

Statistical Analysis
The missing values of climate data and vegetation flux observation data were linearly or non-linearly interpolated using their seasonal dynamics. Based on the ground phenological observations in Haibei and Inner Mongolia, the accuracy of remote sensing simulated phenological results, including MODIS product (m), threshold algorithm derived from NDVI and EVI time series (tn and te), and relative change rate algorithm derived from NDVI and EVI time series (rn and re), were investigated by analysis of variance (ANOVA) and multiple comparisons (Tukey's HSD) (α = 0.05). Linear regression and two statistical parameters, root mean square error (RMSE) and relative prediction error (RPE), were used to quantify the deviation of remote sensing simulation phenology products: where xi is the ground phenological observation, yi is the remote sensing simulated phenological result depending on different methods, and and are the averages of corresponding data, respectively. The n is the number of samples. Root mean square error (RMSE) is used to measure the bias from the simulated data compared to field measurements. The relative predictive error (RPE) provides the direct effectiveness (underestimation as negative, or overestimation as positive) in predicted values compared to measured values [53]. All statistical and modeling procedures were performed in the R statistical computing package (version 3.5.1).

Vegetation Dynamic
The seasonal variation of GPP and vegetation index (NDVI and EVI) in HB and NM grassland ecosystems showed a single-peak curve change pattern (Figures 3a and 4). The maximum generally appeared at the end of July and early August (DOY200-220), and the productivity and vegetation indexes of HB were significantly higher than those of NM (Figures 3a and 4). The interannual dynamics of productivity in HB and NM both showed a significant upward trend during the study period (HB:    (Figures 3a and 4), but the remote sensing estimation of interannual vegetation dynamics was relatively rough (Figures 3b and 4). Specifically, from 2003 to 2011 in HB, the annual integrated vegetation index did not show a significant increase (NDVI: R 2 < 0.1, EVI: R 2 = 0.16, p > 0.05) despite a significant increase of GPP (Figure 3b). In addition, the divergence in interannual dynamics between vegetation index and GPP was more significant in NM. For example, the vegetation index showed a significant downward trend (NDVI: R 2 = 0.9, EVI :R 2 = 0.86, p < 0.05) but GPP gradually increased from 2004 to 2011 (Figure 4).

Vegetation Phenology Dynamic
In this study, a total of six sets of phenological data in HB and NM were obtained using different data and methods (Table 2), which were the phenology measurements on the ground, MODIS products (m), phenological products derived from NDVI and EVI using the threshold method (tn and te) (Figure 2a,c,e,g), and phenological products calculated from NDVI and EVI using the relative rate of change method (rn and re) ( Figure  2b,d,f,h, and Table 1). * Field, m, tn, te, rn, and re represent phenology measurements on the ground, MODIS phenology product (m), phenological products derived from NDVI and EVI using the threshold method (tn and te), and the relative rate of change method (rn and re), respectively. ** Different letters in the same column indicate significant differences in phenological estimation by different methods (p < 0.05) in HB and NM, respectively. *** The negative value of RPE means that the remote sensing phenological estimation method overestimates the measured values.Different phenological algorithms based on remote sensing indexes significantly overestimated the SOS of HB and NM (p < 0.05) ( Table 2), indicating that the phenological simulation based on the remote sensing index of these two sites was not sensitive to SOS. The SOS estimations based on the threshold algorithm had the greatest magnitude of overestimation, especially in HB, where the overestimation reached more than half of the actual value (RPE > 50%, RMSE > 50). Although the phenology estimation calculated from the relative change rate algorithm was relatively improved, there was still a nonnegligible overestimation (RPE < 40%, RMSE < 50). For EOS estimates, the two sites showed different characteristics. For example, the EOS estimations of HB were all overestimated, but the overestimation was obviously smaller than the estimation deviation of SOS (RPE < 10%, RMSE < 20). On the contrary, the EOS simulation based on remote sensing index could provide an unbiased estimate of EOS in NM grassland (p > 0.05, RPE < 5%), especially using the relative change rate algorithm based on EVI data (NM_EOS_re), which could basically accurately estimate the phenological period of NM grassland. (RPE = 0.6%, RMSE = 6.92) ( Table 2).
The interannual variation of the actual observed phenological periods (SOS and EOS) in HB and NM could be well simulated by the five different phenological algorithms based on the remote sensing index introduced in this study (b > 0, Figure 4). However, the accuracy of the trend estimation was different. Specifically, in HB grassland, based on NDVI data, the use of threshold and relative rate of change algorithms to simulate the trend of SOS (HB_SOS_tn and HB_SOS_rn) was more accurate (R 2 > 0.5, p < 0.05), while other SOS simulation methods were not significant (p > 0.05) (Figure 5a). In contrast, HB_EOS_tn could not effectively assess the interannual variability of EOS in HB grassland (R 2 = 0.33, p > 0.1) (Figure 5b). In Inner Mongolia grassland, five different phenological algorithms based on remote sensing index could effectively simulate the interannual variation trend of SOS (R 2 > 0.6, p < 0.05) (Figure 5c), but in comparison, only the EOS estimation using the relative change algorithm based on EVI data (NM_EOS_re) could effectively assess the interannual variability of EOS in NM (R 2 = 0.7, p = 0.01) (Figure 5d).  (Table 2).

Carbon Phenology Dynamics
HB_SOS_rn and HB_EOS_rn provided more ideal performance in the estimations of actual phenological values in HB, while NM_SOS_re and NM_EOS_re could more effectively simulate the phenological period of NM (Figure 5a,b). Therefore, we used the actual phenological values in HB and NM to linearly calibrate the value of satellite-based phenology estimations (HB: HB_SOS_rn and HB_EOS_rn (R 2 > 0.64, p < 0.05), and NM: NM_SOS_re and NM_EOS_re (R 2 = 0.75, p < 0.05), respectively ( Figure 6). In addition, the linearly calibrated satellite-based phenology estimations were considered as the correct phenological values in the years when measured phenology values were missing in HB (HB_SOS_c and HB_EOS_c) (Figure 7a) and in NM (NM_SOS_c and NM_EOS_c) ( Figure  7b).  As GPP in both two sites in the non-growing season was close to 0, according to the threshold method, the ratio of GPP in the time of SOS (or EOS) and the maximum value of GPP within the year (GPP_max), namely GPP_SOS (or GPP_EOS)/GPP_max, was the threshold of carbon phenology. Therefore, based on similar Savitzky-Golay filtering, daily GPP measurements were smoothed at seasonal scales from 2003 to 2011 in HB (GPP_HBs) and in NM (GPP_NMs) (Figure 7c,d). The maximum value of GPP within the year (GPP_max) was obtained from the smoothed GPP temporal series and combined with the phenological values determined from the actual phenological observations (or the linearly calibrated satellite-based phenology estimations) (Figure 7a,d) to estimate the threshold of carbon phenology in HB and NM, respectively.
The carbon phenology thresholds of grasslands both in HB and NM showed obvious interannual variations (Figure 8). Specifically, during the observation period, the carbon phenology (SOS and EOS) thresholds of HB grassland (HB_SOSt and HB_EOSt) had a significant upward trend, with the multiyear average values being 0.14 and 0.29, respectively. HB_SOSt reached its maximum value in 2010 (0.25). HB_SOSt was relatively smaller in the first three years at the beginning of the observation period, with an average value of 0.02. In contrast, the thresholds of Inner Mongolia grassland (NM_SOSt and NM_EOSt) also had a certain upward trend, but it was not significant (p > 0.05), with the multiyear average values being 0.17 and 0.2, respectively. In addition, the interannual variability of the carbon phenology thresholds in NM was larger than that of HB grassland, with the exception of 2009 when both NM_SOSt and NM_EOSt exceeded 0.5. For the rest of the years, NM_SOSt and NM_EOSt were less than 0.3 ( Figure 8).

Satellite-based Vegetation Dynamics
Satellite-based vegetation indexes could match the seasonal dynamics of vegetation growth well, but they were relatively coarse in their estimation of interannual vegetation dynamics (Figures 2 and 3). Most studies also found that satellite-based vegetation indexes and GPP models could capture the intraannual vegetation dynamics well [54][55][56][57][58]. However, assessing interannual vegetation dynamics using satellite-based estimations is full of challenges because of complex interactions between biotic and environmental conditions [56,59,60]. Plant phenology is mainly dependent on intraannual recurring plant growth and reproductive phenomena [61]. Therefore, interannual mismatch between satellite-based vegetation indexes and field observations does not interfere with phenological extraction from satellite-based vegetation dynamics.

Satellite-based Vegetation Phenology Estimations
Satellite-based vegetation phenology estimations (all five of the algorithms in this study) strongly coupled with interannual variations of the actual observed phenological period ( Figure 5) but significantly overestimated the SOS in HB and NM (>20% of RPE, Table 2). In other words, the plants had actually entered their growing season, but the satellite-based vegetation phenology estimations did not detect this. Many studies also indicated consistent trends in vegetation phenology between satellite-based phenology estimations and field or near-surface remote observations [22,24,31,32,40,62]. Confounding interferences from the background reflectance of soil, snow, leaf litter, and shadows may lead to these overestimations in satellite-based vegetation phenology [30]. For example, mounds of leaf litter and withered grass in these cold temperate grasslands in the early growing season masked the fact that the plants were turning green [63]. This interference from leaf litter weakened as the vegetation grew taller; thus, satellite-based EOS estimations had a better performance at both sites than SOS estimations (Table 2). Additionally, snow cover and snowmelt remarkably influenced the effect of satellite-based phenology algorithms due to the similar sensitivity of remote-based vegetable indexes to plant greening and snowmelt [19,64,65]. Therefore, field observations and the use of nearsurface remote sensing are still necessary to evaluate and refine satellite-based phenology estimates.

Vegetation Carbon Phenology and Thresholds
Plant carbon phenology is strongly linked to but not completely equivalent to the natural vegetation phenology period [41][42][43]66,67]. For example, satellite-based SOS estimates from MODIS were highly correlated with the onset of photosynthesis derived from flux measurements despite a later bias for coniferous sites [43]. However, the accuracy in describing annual patterns of flux phenology was generally coarse using satellite-based phenology estimates from different data resources [42]. Thus, most studies found that the trend of vegetation carbon phenology and plant phenology was consistent, but those of absolute dates were not well matched [3,26,42,43]. One reason was due to the uncertainty of satellite-based vegetation phenology estimations and the data quality of EC-based observations [64]. Another possible reason was that the empiric prethreshold settings were not fixed on the space-time scale [68]. Actually, our results showed that there was considerable interannual variation in carbon phenology thresholds in two cold temperate grasslands in China ( Figure 8).
Statistically, the EC-based maximum GPP observations (GPPmax) of the year had a significant correlation with the thresholds of carbon phenology periods (SOSt and EOSt) at the HB site (Figure 9a), which provided a simple and valid way to calculate SOSt and EOSt at a specific site and time. Further, the interannual variations of SOSt and EOSt were weakly correlated (Figure 9b), which was largely due to the influence of spring vegetation green-up dates on EOS [1]. The relationship between SOSt and EOSt also implied an alternative way to analyze EOSt variations based on SOSt variations due to a limit of predictive strength in EOS estimation in most phenology studies [1,26,33]. Therefore, this study provided bridges between satellite-based vegetation phenology and EC-based carbon phenology. However, many uncertainties that we did not consider in this study, i.e., agreement of EC-based GPP products with a large footprint and local-scale ground measurements, and the phenology match of specific species and grassland ecosystems across different ecosystems, may have had confounding effects on vegetation carbon phenology and thresholds estimations. Thus, spatial extrapolation to different vegetation types or regions should be based on appropriate validation. Figure 9. The determination of the thresholds of carbon phenology periods from EC-based maximum GPP observations (GPPmax) in HB (a) and the relationship between two thresholds of carbon phenology periods (SOSt and EOSt) in HB and NM sites (b). The solid line represents a significant regression relationship (p < 0.05).

Conclusions
This study investigated vegetation phenology in two typical cold temperate grasslands (HB and NM) in China from 2001 to 2017. Compared to ground vegetation phenology observations, we analyzed the performance of MODIS phenology products and two remote sensing-based vegetation phenology algorithms from NDVI and EVI time series. The vegetation phenology estimation from the optimal algorithm was used to calculate the thresholds of carbon phenology periods (SOSt and EOSt) from EC-based GPP observations in each site. Results showed single-peak seasonal variations and significant upward interannual trends of local vegetation growth during the study period in both sites (HB: R2 = 0.79, NM: R2 = 0.57, p < 0.05). Satellite-based phenology estimations (all five algorithms in this study) strongly coupled with temporal variation of the actual observed phenological period but significantly overestimated the SOS, predicting it to be over 21 days later than the field data. The carbon phenology thresholds of HB and NM grasslands showed obvious interannual variations. Specifically, the carbon phenology (SOS and EOS) thresholds of HB grassland (HB_SOSt and HB_EOSt) had a significant upward trend, with the multiyear average values being 0.14 and 0.29, respectively. In contrast, the thresholds of NM grasslands (NM_SOSt and NM_EOSt) also have a certain upward trend but not significant (p>0.05), with the multiyear average values being 0.17 and 0.2, respectively. Our study suggested the thresholds of carbon phenology periods (SOSt and EOSt) could be simply estimated based on their significant relationship with the EC-based maximum GPP observations (GPPmax) in a specific site and time. However, spatial extrapolation to different vegetation types or regions should be based on appropriate validation.