Spatial Difference between Temperature and Snowfall Driven Spring Phenology of Alpine Grassland Land Surface Based on Process-Based Modeling on the Qinghai–Tibet Plateau

: As a sensitive indicator for climate change, the spring phenology of alpine grassland on the Qinghai–Tibet Plateau (QTP) has received extensive concern over past decade. It has been demonstrated that temperature and precipitation/snowfall play an important role in driving the green-up in alpine grassland. However, the spatial differences in the temperature and snowfall driven mechanism of alpine grassland green-up onset are still not clear. This manuscript establishes a set of process-based models to investigate the climate variables driving spring phenology and their spatial differences. Speciﬁcally, using 500 m three-day composite MODIS NDVI datasets from 2000 to 2015, we ﬁrst estimated the land surface green-up onset (LSGO) of alpine grassland in the QTP. Further, combining with daily air temperature and precipitation datasets from 2000 to 2015, we built up process-based models for LSGO in 86 meteorological stations in the QTP. The optimum models of the stations separating climate drivers spatially suggest that LSGO in grassland is: (1) controlled by temperature in the north, west and south of the QTP, where the precipitation during late winter and spring is less than 20 mm; (2) driven by the combination of temperature and precipitation in the middle, east and southwest regions with higher precipitation and (3) more likely controlled by both temperature and precipitation in snowfall dominant regions, since the snow-melting process has negative effects on the air temperature. The result dictates that snowfall and rainfall should be concerned separately in the improvement of the spring phenology model of the alpine grassland ecosystem.


Introduction
The Earth has experienced significant warming since the industrial revolution (IPCC, 2013), and the warming rate in the Qinghai-Tibet Plateau (QTP) over the past five decades is twice the global average [1]. Rapid warming has significantly affected the biogeochemical cycles [1] and the ecological environment [2] in QTP grassland, which covers 66% of the area. In particular, the warming-deduced variation in vegetation spring phenology directly affects the grassland's net primary productivity and herbage production, which has a strong impact on the development of the regional economy. As a result, the date of green-up onset, which represents the start of the growing season in grasslands on the QTP, has recently been extensively investigated [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. The majority of previous research is based on statistical methods to reveal the relationship between the spring phenology and climate variables [3,6,7,10,12,15,16,19,20], however, process models and process simulations of green-up onset are not well applied, and the climate control mechanism of grassland green-up onset and the spatial difference are not well investigated [11,21]. It has been recognized that a process-based model could effectively and accurately simulate the plant development process and the phenological phase based on various climate variables, such as temperature, precipitation and photoperiod [22]. Therefore, it is important to develop and apply the process-based spring phenology model for scientifically understanding the influence mechanism of climate change on the green-up date of grasslands and analyzing the spatial difference on the QTP for dealing with climate change, ecological protection and animal husbandry development.
Vegetation spring phenology has been commonly modeled by associating ground based phenology observations for individual species with meteorological data [22].This is also true on the QTP [11]. Although temperature is regarded as the main driving factor for the simulation of woody plant spring phenology [23,24], temperature, precipitation and their coupling effect are usually integrated together in the construction of a spring phenology model in grasslands [11,21,25,26]. Specifically, grassland green-up has been linked to air temperature, precipitation and their coupling effect as observed in meteorological stations [3][4][5]10,27], suggesting that precipitation or snowfall have much stronger controls on grassland green-up [7,11,28]. Moreover, a process-based phenology model, which was established based on multi-stations, multi-species and long term series of herb green-up dates from ground observations in the east and south of the QTP, indicates that late winter and spring snowfall is an important factor that affects the herb green-up [11]. However, the ground phenology observations are very limited and are mainly distributed in the margin of the QTP [11,29,30], so that the resultant phenological models only represent the corresponding phenological sites, but they are generally unable to explain detailed spatial patterns of climate controls on alpine grassland spring phenology on a large scale.
Specifically, varieties of models have been proposed to simulate the vegetation spring phenology [31], such as the GDD (growing degree days) model [32][33][34], sequential model, parallel model, alternating model and the unified model [23], with temperature as the only driving factor. However, these models are mainly developed for the simulation of spring phenology (leaf unfolding/leaf-out/budburst) for woody plants. For herbaceous plants with different physiological characteristics, because of the existing precipitation's control over the spring phenology, precipitation was embedded in the above-mentioned models to improve the simulation of the green-up date [11,25]. However, either for woody plants or herbaceous plants, the models were designed based on the physiological processes of individual plants to simulate plant phenology at the species-level, which is dependent on limited plant individual phenology observation and could not be generalized to the land surface in a large area.
In contrast, land surface phenology derived from satellite data provides exhaustive spatial distributions, which can be used to establish much more robust phenology models to illustrate climate impacts on phenology variation on the QTP [21,35]. However, current detections of land surface phenology vary greatly, with different methods of processing the noisy time series of satellite observations [36,37]. Most of the phenology detections contain large uncertainties compared with ground observations [37]. Thus, it is necessary to evaluate the process-based simulation for the land surface green-up onset (LSGO) based on in situ observations.
In this study, we aim to investigate the following scientific questions: (1) can LSGO in grassland be used for the establishment of process-based models for exploring the spatial variation in climate drivers across QTP? (2) Can process-based models reveal the spatial difference of precipitation/snowfall control on alpine grassland spring phenology?

Study Area and Biomes
As the highest plateau in the world, the QTP is called 'the roof of the world', with a mean elevation exceeding 4000 m. Due to its special alpine environment, the QTP is the largest glacier store region, following the Arctic and Antarctica, and is called 'the third pole' of the Earth. The QTP is mainly distributed in China, including Tibet Autonomous Region, western Sichuan Province, part of Yunnan Province, Qinghai Province, southern Xinjiang Uygur Autonomous Region and part of Gansu province. The QTP climate is affected by south Asia monsoon in summer and prevailing westerlies in winter. From northwest to southeast, perennial mean air temperature ranges from −5 to 15.5 • C and multi-year average precipitation has been 16-1764 mm during the past three decades. Unique vegetation types have evolved on the QTP due to its special physic-geographical environment (Figure 1), which mainly consists of herbaceous plants that are steppe, meadow and alpine vegetation over 66% of the area [38]. Specifically, steppe is mainly dominated by Carex grass high-cold steppe, and the other types (which are only distributed in temperate regions on the northern QTP), include temperate tufted grass steppe, temperate tufted low grass nano-semi-shrub desert steppe and temperate grass forb meadow steppe [38]. Meadow is mainly dominated by Kobresia forb high-cold meadow, the other types including Carex grass and forb swamp meadow, grass and forb halophytic meadow and grass and forb meadow [38]. Alpine vegetation mainly consists of the composite family including alpine tundra and sparse vegetation and alpine cushion vegetation [38]. The other vegetation types include broad-leaved and coniferous forests in the south and southeast and shrub and cultivated vegetation in the east and margin of the east, but the areas are smaller [38].
As the highest plateau in the world, the QTP is called 'the roof of the world', with a mean elevation exceeding 4000 m. Due to its special alpine environment, the QTP is the largest glacier store region, following the Arctic and Antarctica, and is called 'the third pole' of the Earth. The QTP is mainly distributed in China, including Tibet Autonomous Region, western Sichuan Province, part of Yunnan Province, Qinghai Province, southern Xinjiang Uygur Autonomous Region and part of Gansu province. The QTP climate is affected by south Asia monsoon in summer and prevailing westerlies in winter. From northwest to southeast, perennial mean air temperature ranges from -5 to 15.5 ℃ and multiyear average precipitation has been 16-1764 mm during the past three decades. Unique vegetation types have evolved on the QTP due to its special physic-geographical environment (Figure 1), which mainly consists of herbaceous plants that are steppe, meadow and alpine vegetation over 66% of the area [38]. Specifically, steppe is mainly dominated by Carex grass high-cold steppe, and the other types (which are only distributed in temperate regions on the northern QTP), include temperate tufted grass steppe, temperate tufted low grass nano-semi-shrub desert steppe and temperate grass forb meadow steppe [38]. Meadow is mainly dominated by Kobresia forb high-cold meadow, the other types including Carex grass and forb swamp meadow, grass and forb halophytic meadow and grass and forb meadow [38]. Alpine vegetation mainly consists of the composite family including alpine tundra and sparse vegetation and alpine cushion vegetation [38]. The other vegetation types include broad-leaved and coniferous forests in the south and southeast and shrub and cultivated vegetation in the east and margin of the east, but the areas are smaller [38].

Ground Observations of Grass Green-Up and Meteorological Datasets
The grass green-up dates observed in fields are obtained from the China Meteorological Administration, which consists of 19 herb species over 6 phenological ground stations from 2000 to 2012. There are two to eight herb species being observed in each station (Table 1). The green-up date is observed by professional staff based on the same plant observation method with a united observation specification [39]. Specifically, the typical grassland in each station is selected, and then a 0.5 m × 0.5 m grassland quadrat is randomly extracted. The growing stages of 10 individuals for each species within the quadrat are observed every 2 days. The green-up date of the herbs is defined as the date when the young leaves in 50% of individuals grow higher than 1 cm over the ground. Usually, the

Ground Observations of Grass Green-Up and Meteorological Datasets
The grass green-up dates observed in fields are obtained from the China Meteorological Administration, which consists of 19 herb species over 6 phenological ground stations from 2000 to 2012. There are two to eight herb species being observed in each station ( Table 1). The green-up date is observed by professional staff based on the same plant observation method with a united observation specification [39]. Specifically, the typical grassland in each station is selected, and then a 0.5 m × 0.5 m grassland quadrat is randomly extracted. The growing stages of 10 individuals for each species within the quadrat are observed every 2 days. The green-up date of the herbs is defined as the date when the young leaves in 50% of individuals grow higher than 1 cm over the ground. Usually, the observations are focused on the dominant species, such as poaceae, leguminosae and cyperaceae, and accompanying species, such as compositae, chenopodiaceae, liliaceous, rosaceae, umbrella and labiatae ( Table 1). The mean green-up dates from multiple species in each station are used to establish the process-based models to evaluate the accuracy of land surface green-up onset (LSGO) (cf. Section 2.4). The meteorological datasets are derived from China's meteorological data sharing service system in the China Meteorological Administration (http://cdc.nmic.cn/home.do (accessed on 22 December 2021)). The data include daily mean air temperature and daily precipitation from 2000 to 2015 in 86 stations ( Figure 1). The stations are located at the elevation range from 1813.9 to 4900 m. Since the snowfall and rainfall data in late winter and spring (from 1 January to the multi-year mean green-up date) cannot be directly acquired for the meteorological stations, the snowfall data are defined as the accumulation of precipitation with daily mean air temperature below 0 • C [40]. The total rainfall is the accumulation of precipitation with daily mean air temperature above or equal to 0 • C.

Land Surface Green-Up Onset Dataset and Evaluation
Land surface spring phenology from 2000 to 2015 was estimated based on the Moderate Resolution Imaging SpectroRadiometer (MODIS) Nadir Bidirectional Reflectance Distribution Function (BRDF)-Adjusted Reflectance (NBAR) dataset (MCD43A4 Version 6) and Albedo Quality Product (MCD43A2) with a spatial resolution of 500 m and a temporal resolution of one day. Both MCD43A4 and MCD43A2 products were downloaded from the website (http://modis.gsfc.nasa.gov/ (accessed on 22 December 2021)) of the National Aeronautics and Space Administration (NASA). Daily normalized differential vegetation index (NDVI) was calculated from MCD43A4.
Three-day composite NDVI was computed first using maximum principle. Then, NDVI time series were fitted using the vegetation growing curve described in a logistic model, and curvature change rate was calculated to identify the phenological transition dates [41,42]. Specifically, a Hybrid Piecewise Logistic Model (HPLM) algorithm was first used to describe the temporal NDVI trajectory, optimizing the parameters of the logistic equation. Then, the curvature of the equation was computed. Finally, curvature variation (first derivative of curvature) was computed and the independent variable (day of year) which corresponds to the first extreme value of the dependent variable (first derivative of curvature) was defined as the land surface green-up onset (LSGO) (Figure 2). Details of the method are given in Zhang's paper [43,44]. It should be note that the first growing cycle of grassland is only addressed in this study if there are two or more growing cycles. The land surface green-up onset dataset was evaluated using the multi-species mean green-up date in the six stations by mean absolute error (MAE) (Equation (1)), root mean square error (RMSE) (Equation (2)) and Pearson correlation coefficient (r) (Equation (3)), which are described as the followings.
where i is the sequence number of the year, n is the total number of years, x i is the mean green-up date of ground observations from multi species in the i-th year, y i is LSGO in the corresponding pixel, x is the mean of multi-year ground observations in multi-species and y is the mean of multi-year LSGO in one or multi-pixels.
In the evaluation, considering the uncertainties of remote sensing surface reflectance caused by the outliers and the spatial heterogeneity for the land surface phenology, field observations in each station are compared with LSGO within a set of pixel-window sizes to identify the optimum window size according to MAE, RMSE and r. The pixel-window size varies from one to twenty-one 500 m-pixels with the field station located in the center [45]. We then calculated the medium green-up onset of pixels within 9.5 km 2 centered with the meteorological station in each window size to present the green-up onset of the alpine grassland of the meteorological station, which could also ensure the meteorological data are representative within this extent (usually less than 10 km 2 [45]).
(first derivative of curvature) was computed and the independent variable (day of year) which corresponds to the first extreme value of the dependent variable (first derivative of curvature) was defined as the land surface green-up onset (LSGO) (Figure 2). Details of the method are given in Zhang's paper [43,44]. It should be note that the first growing cycle of grassland is only addressed in this study if there are two or more growing cycles. The land surface green-up onset dataset was evaluated using the multi-species mean green-up date in the six stations by mean absolute error (MAE) (Equation (1)), root mean square error (RMSE) (Equation (2)) and Pearson correlation coefficient (r) (Equation (3)), which are described as the followings.
where i is the sequence number of the year, n is the total number of years, is the mean green-up date of ground observations from multi species in the i-th year, is LSGO in the corresponding pixel, ̅ is the mean of multi-year ground observations in multi-species and is the mean of multi-year LSGO in one or multi-pixels.
In the evaluation, considering the uncertainties of remote sensing surface reflectance caused by the outliers and the spatial heterogeneity for the land surface phenology, field observations in each station are compared with LSGO within a set of pixel-window sizes to identify the optimum window size according to MAE, RMSE and r. The pixel-window size varies from one to twenty-one 500 m-pixels with the field station located in the center [45]. We then calculated the medium green-up onset of pixels within 9.5 km 2 centered with the meteorological station in each window size to present the green-up onset of the alpine grassland of the meteorological station, which could also ensure the meteorological data are representative within this extent (usually less than 10 km 2 [45]). .

Establishment of Process-Based Model for Green-Up Dates
Process-based models for vegetation green-up date are established using daily mean air temperature and precipitation. The models are mainly divided into three groups [11,25], which are the traditional thermal time model, air temperature-precipitation parallel model and air temperature-precipitation sequential model. The traditional thermal time model as-sumes that vegetation begins to green-up when the accumulated temperature R T (T t ) reaches to a critical value F * , which is described as (Equation (4)): where t 0 is the initial day, which is usually defined as 1 January, y is the date of green-up onset and S T is the state of forcing temperature. R T (T t ) is defined as Equations (5) and (6). T b is the degree-day base temperature indicating the linear response of green-up to the accumulation of temperature, a and b are the parameters and T t is the daily mean air temperature. Equation (5) is commonly used in previous studies [11,25]. However, the response of plant to temperature is usually nonlinear [46]. Thus, in our research, Equation (6) with logistic response function was used in the construction of the spring phenology model. The air temperature-precipitation parallel model assumes that vegetation begins to green-up when the accumulated temperature reaches to a critical value F * (Equation (4)) and the state of forcing in precipitation (S p ) reaches a critical value P * from a date of t 0 (Equations (7) and (8)): where S p is the state of forcing in precipitation, P t is the daily precipitation (mm) and R p (P t ) is the state of precipitation as a daily sum of the rate of precipitation. The air temperature-precipitation sequential model assumes that the effect of air temperature on green-up (Equation (4)) is enforced after the state of forcing in precipitation reaches a critical value P * at the date of t1 (Equation (9)), which means that, after this date, the forcing effect of air temperature on herb green-up will work until the date of green-up onset (y): Four parameters of a, b, F * and P * should be determined for the air temperatureprecipitation parallel and sequential models. The parameters of the three models are fitted using the simulated annealing algorithm of Metropolis [47] and are constrained by the lowest RMSE (root mean square error (Equation (10)) between simulated and observed green-up dates: where obs i is the observed green-up date in year i, sim i is the simulated green-up date in year i and n is the number of the year. RMSE is a common indicator for the model evaluation. Usually, lower RMSE means better simulation accuracy. The aforementioned three process-based models are established for ground observed grass green-up dates in each of six agrometeorological stations and for LSGO in each of eight-six meteorological stations. An optimum model for each station is chosen from the three process-based models to identify the factors driving green-up dates. The optimum model is selected based on the lowest value of Akaike information criterion (AIC) [11,25] (Equation (11)), which rates models in terms of parsimony and efficiency: where obs i is the observed green-up date in year i, sim i is the simulated green-up date in year i, n is the number of the year and k is the number of the parameters. It is a commonly used index for the evaluation of the process-based phenology models [31,48], which could effectively avoid overfitting. It should be noted that the 16-year time series of LSGO in our research is enough for the phenology process-based simulation (usually more than 10 years) [11,25], so the data after 2015 are not included in our research.

Land Surface Green-Up Onset of Alpine Grassland and Evaluation
The multi-year mean LSGO in grasslands varies from March to July (Figure 3). It shows that the green-up onset of alpine grassland tends to be later from the east to the west of the QTP. This spatial pattern is same as the results in previous studies [6,49], which demonstrates the high reliability of the MODIS green-up onset dataset. The MAE and RMSE of the 6 stations in different pixel windows are less than 20 days and the r of 4 stations shows significant positive correlations ( Figure A1), which further demonstrated high reliability of the dataset. eight-six meteorological stations. An optimum model for each station is chosen from the three process-based models to identify the factors driving green-up dates. The optimum model is selected based on the lowest value of Akaike information criterion (AIC) [11,25] (Equation (11)), which rates models in terms of parsimony and efficiency: = × ∑ + 2 + 1 (11) where is the observed green-up date in year i, is the simulated green-up date in year i, n is the number of the year and k is the number of the parameters. It is a commonly used index for the evaluation of the process-based phenology models [31,48], which could effectively avoid overfitting. It should be noted that the 16-year time series of LSGO in our research is enough for the phenology process-based simulation (usually more than 10 years) [11,25], so the data after 2015 are not included in our research.

Land Surface Green-Up Onset of Alpine Grassland and Evaluation
The multi-year mean LSGO in grasslands varies from March to July (Figure 3). It shows that the green-up onset of alpine grassland tends to be later from the east to the west of the QTP. This spatial pattern is same as the results in previous studies [6,49], which demonstrates the high reliability of the MODIS green-up onset dataset. The MAE and RMSE of the 6 stations in different pixel windows are less than 20 days and the r of 4 stations shows significant positive correlations ( Figure A1), which further demonstrated high reliability of the dataset.  Table 2 presents the optimum models that are selected for the ground green-up date and the corresponding LSGO. In six agrometeorological stations (both phenological and meteorological data are available), the optimum model in each station is the best one to simulate green-up dates among the traditional thermal time model, air temperature-precipitation parallel model and air temperature-precipitation sequential model. The results indicate that the same type of optimum model is selected for both ground and LSGO  Table 2 presents the optimum models that are selected for the ground green-up date and the corresponding LSGO. In six agrometeorological stations (both phenological and meteorological data are available), the optimum model in each station is the best one to simulate green-up dates among the traditional thermal time model, air temperatureprecipitation parallel model and air temperature-precipitation sequential model. The results indicate that the same type of optimum model is selected for both ground and LSGO green-up dates in two agrometeorological stations (Henan and Gade). However, all the optimal models for both ground observations and LSGO detections show that the green-up date is driven by the combination of air temperature and precipitation together (either air temperature-precipitation parallel or sequential model), except for the field green-up date in Xinghai, where the traditional thermal time model works better. Moreover, the RMSEs are very close to each other for ground observations and the LSGO based optimum model. The results overall suggest that the LSGO model is capable of simulating the grass green-up date with a similar accuracy as the ground based model. F* and P* are critical values of state of forcing from quiescence to green-up in temperature and precipitation, respectively, M t is traditional thermal time model, M t-p1 is air temperature-precipitation parallel model, M t-p2 is air temperature-precipitation sequential model, multi-species is the multi-species mean green-up date of herbs in each station and multi-pixel is the median SGO within 9.5 km square centered the stations. * p < 0.05, ** p < 0.01, *** p < 0.001. Figure 4a presents the optimal process-based models in each of 86 meteorological stations established based on medium LSGO in a window size of 9.5 km and daily mean air temperature and precipitation. The optimum models present clear spatial clusters geographically (Figure 4a). The thermal time model appears in 23 stations (27%) and is mainly distributed in the north, west and south of the QTP. This indicates that temperature drives green-up of grassland for the stations. In contrast, the temperature-precipitation parallel model (41% of stations) and temperature-precipitation sequential model (32% of stations) are spatially mixed and mainly distributed in the middle, east and southwest regions. This suggests that the green-up onset of grassland is mainly driven by the combination of air temperature and precipitation.

The Spatial Pattern of Precipitation/Snowfall in the Stations
The annual mean preseason precipitation in each station spreads from 4.8 to 365 mm. On the whole, the precipitation in the south and southeast stations is higher than other stations, while the spatial difference is not obvious (Figure 4b). However, the annual mean preseason snowfall of the 86 stations shows obvious spatial difference (Figure 4c). Specifically, the stations with snowfall less than 20 mm are mainly distributed in the north, west and south QTP, while the stations with snowfall higher than 20 mm are mainly distributed in middle, east and southwest edge QTP. Furthermore, the stations with a proportion of snowfall (ration of snowfall and precipitation) higher than 50% were mainly distributed in the middle, east and southwest edge of the QTP (Figure 4b).

The Spatial Pattern of Precipitation/Snowfall in the Stations
The annual mean preseason precipitation in each station spreads from 4.8 to 365 mm. On the whole, the precipitation in the south and southeast stations is higher than other stations, while the spatial difference is not obvious (Figure 4b). However, the annual mean preseason snowfall of the 86 stations shows obvious spatial difference (Figure 4c). Specifically, the stations with snowfall less than 20 mm are mainly distributed in the north, west and south QTP, while the stations with snowfall higher than 20 mm are mainly distributed in middle, east and southwest edge QTP. Furthermore, the stations with a proportion of snowfall (ration of snowfall and precipitation) higher than 50% were mainly distributed in the middle, east and southwest edge of the QTP (Figure 4b).

Reliability of MODIS Land Surface Green-Up Onset for Construction Process-Based Model on the QTP
The process-based model for LSGO is expected to provide monitoring and prediction of phenology in vegetation communities across large scales. This could overcome the drawbacks in commonly used process-based phenology models established using extremely limited in situ observations. The process-based model using LSGO has shown advantages over a large region [35,50], such as the model for the woody leaf unfolding date in the deciduous broadleaf forest region of China [50].
However, although the MODIS LSGO has been used for the process simulation of alpine grassland green-up on the QTP [21], it is still full of uncertainty [36,37] and shows

Reliability of MODIS Land Surface Green-Up Onset for Construction Process-Based Model on the QTP
The process-based model for LSGO is expected to provide monitoring and prediction of phenology in vegetation communities across large scales. This could overcome the drawbacks in commonly used process-based phenology models established using extremely limited in situ observations. The process-based model using LSGO has shown advantages over a large region [35,50], such as the model for the woody leaf unfolding date in the deciduous broadleaf forest region of China [50].
However, although the MODIS LSGO has been used for the process simulation of alpine grassland green-up on the QTP [21], it is still full of uncertainty [36,37] and shows low predicted ability [21], which needs further evaluation based on in situ observations. Based on ground herb green-up date observation in six alpine grassland stations at the middle and east QTP, we evaluated the MODIS LSGO. Since the ground herb observation (tens of meters) does not spatially match with the MODIS pixel scale (hundreds of meters, community scale), it is hard to directly compare their difference [51]. However, we found that the MODIS LSGO could capture the green-up date variation in most of the stations. The simulation results from using both ground herb green-up date and LSGO at the same stations, such as lower RMSE and higher R 2 , supported our hypothesis that LSGO could perform well for the construction of process-based spring phenology. Specifically, the LSGO model simulation is similar to the multi-species mean green-up date observed from 33% of ground stations. However, the driving factors derived from the LSGO optimum model simulation are consistent with those from simulation of optimum models of the ground observations in 83% of stations. This suggests that the process-based optimum LSGO model could reflect well for the green-up onset climate driving factors in grassland, which is potentially capable of revealing the alpine grassland spring phenology climatic driving mechanism across the QTP.

Spatial Difference of the Optimum Process-Based Model for LSGO
The optimum LSGO model distinguishes the different climate drivers of grass greenup dates across 86 meteorological stations. The thermal time models and temperatureprecipitation models reveal that both precipitation and temperature together trigger the grass green-up date in the middle, east and southwest regions, but thermal time alone drives the grass green-up timing in north, west and south areas (Figure 4a). This spatial difference in phenological forces is likely also associated with the amount of snowfall. A previous study based on the process-based model of the in situ herb green-up date shows that the green-up is mainly driven by air temperature in the stations with less than 20 mm snowfall, while it is mainly controlled by the coupling process between air temperature and precipitation in other regions [11]. Similarly, the optimum LSGO model in this study reveals that LSGO is manly driven by air temperature in the areas where precipitation in late winter and spring is less than 20 mm, while the LSGO is manly driven by air temperature and precipitation in other areas where precipitation is relatively high (Figure 4a,b). Since the regions with precipitation of less than 20 mm are mostly covered by less than 20 mm snowfall (Figure 4c,d), this suggests that air temperature is the main driver in triggering grass green-up in the regions with both precipitation and snowfall of less 20 mm. Nevertheless, the snow melting in the region with high snowfall is likely to reduce the air temperature, which could complicate climate controls on grass green-up. In this case, air temperature-precipitation models should be more appropriate for simulating grass green-up timing.

The Role of Snowfall in Driving Alpine Grassland Spring Phenology
Although temperature was regarded as the primary climatic factor controlling the plant spring phenology [52], the in situ species-level observations [28] and manipulative experiment [16] in alpine meadow on the QTP have demonstrated that soil moisture is also important in controlling green-up onset of alpine grassland on the QTP [28,53]. With the decrease in the soil moisture in spring, the green-up of the vegetation community was significantly delayed as the drought was induced [54]. Snowmelt is the only source of soil water in winter [28] in many areas of the QTP. By affecting soil moisture, snowfall could further control the alpine herb green-up. Moreover, based on multi-source remote sensing data streams, it has been reported that the snow cover melt date is positively correlated to the green-up onset in most areas of the QTP [17,27,55]. However, our research confirmed the snow control on alpine grassland on the QTP, but in limited spatial coverage. Except for the effect on soil moisture, the snow also has an indirect impact on the green-up onset of alpine grassland by reducing temperature when it melts [11]. The indirect control of snow on the green-up onset varies with the amount of snowfall in the late winter and early spring. The process-based model simulations in our research demonstrated that the greenup onset at the stations (80.7%) with larger amounts of snowfall (usually more than 20 mm) is mainly controlled by temperature and snowfall (31 stations in total, 25 stations' optimum models are the parallel or sequential model, Figure 4a,c), which support the conjunction effect of snowfall and temperature. However, the dry winter with less snowfall in the alpine steppe distribution areas on the western and northern QTP resulted in temperature controlled green-up. Moreover, the proportion of snowfall in preseason precipitation is also an important factor that should be investigated. A higher proportion of snowfall means stronger impact on the temperature when the snowfall melts, which indicated a stronger control on the alpine grassland green-up onset. In fact, the conjunction effect of snowfall and temperature made the climate control on the alpine grassland more complex, which Remote Sens. 2022, 14, 1273 11 of 14 highlights that there should be different control degrees between preseason snowfall and rainfall on alpine grassland green-up onset, suggesting the snowfall and rainfall should be considered differently in the phenological module in the ecosystem model.

Conclusions
The optimum process-based models of the land surface green-up date are functionally consistent with the models of ground observed grass green-up date in 83% of the stations. This suggests that MODIS land surface green-up date is applicable for the construction of spring phenology process-based models and captures well the inter-annual variation in the ground green-up date in alpine grassland in the QTP. Thus, the process-based model of the land surface green-up date could be applied for monitoring grass growth across the QTP.
The spatial variation in optimum process-based models of land surface green-up date clearly distinguishes the different climate drivers controlling the grass green-up date across the QTP. The thermal time model is dominated in the north, west and south of the QTP, which implies that the green-up timing in grassland is mainly triggered by the spring temperature. In contrast, temperature-precipitation parallel and sequential models are mainly distributed in the middle, east and southwest, where the grassland green-up date is driven by the combination of temperature and precipitation.
Finally, the grassland green-up on the QTP could also be influenced by snowfall. The impact from the snowfall and rainfall to the grassland green-up is different because the snow melting process has negative effects on the air temperature. Future research is needed to separate snowfall and rainfall during late winter and spring in process-based phenological models.  Acknowledgments: The authors thank the Meteorological Information Center of the China Meteorological Administration for providing phenological and meteorological datasets. We also thank NASA for providing free and ready access to the MODIS data.

Conflicts of Interest:
The authors declare no conflict of interest.