Improved Modeling of Gross Primary Productivity of Alpine Grasslands on the Tibetan Plateau Using the Biome-BGC Model

: The ability of process-based biogeochemical models in estimating the gross primary productivity (GPP) of alpine vegetation is largely hampered by the poor representation of phenology and insu ﬃ cient calibration of model parameters. The development of remote sensing technology and the eddy covariance (EC) technique has made it possible to overcome this dilemma. In this study, we have incorporated remotely sensed phenology into the Biome-BGC model and calibrated its parameters to improve the modeling of GPP of alpine grasslands on the Tibetan Plateau (TP). Speciﬁcally, we ﬁrst used the remotely sensed phenology to modify the original meteorological-based phenology module in the Biome-BGC to better prescribe the phenological states within the model. Then, based on the GPP derived from EC measurements, we combined the global sensitivity analysis method and the simulated annealing optimization algorithm to e ﬀ ectively calibrate the ecophysiological parameters of the Biome-BGC model. Finally, we simulated the GPP of alpine grasslands on the TP from 1982 to 2015 based on the Biome-BGC model after a phenology module modiﬁcation and parameter calibration. The results indicate that the improved Biome-BGC model e ﬀ ectively overcomes the limitations of the original Biome-BGC model and is able to reproduce the seasonal dynamics and magnitude of GPP in alpine grasslands. Meanwhile, the simulated results also reveal that the GPP of alpine grasslands on the TP has increased signiﬁcantly from 1982 to 2015 and shows a large spatial heterogeneity, with a mean of 289.8 gC / m 2 / yr or 305.8 TgC / yr. Our study demonstrates that the incorporation of remotely sensed phenology into the Biome-BGC model and the use of EC measurements to calibrate model parameters can e ﬀ ectively overcome the limitations of its application in alpine grassland ecosystems, which is important for detecting trends in vegetation productivity. This approach could also be upscaled to regional and global scales.


Introduction
Gross Primary Productivity (GPP), defined as the total amount of carbon fixed by vegetation through photosynthesis, is one of the most important components of the terrestrial carbon cycle [1,2]. Meanwhile, it is also the largest carbon flux between the terrestrial biosphere and the atmosphere, and drives a variety of ecosystem services such as food, fiber, and wood production [3]. The Tibetan Plateau (TP), with an average elevation exceeding 4000 m, is one of the most sensitive regions to global climate change [4]. Moreover, due to the positive feedbacks associated with the melting of snow and parameter calibration method for complex process-based models; (3) to analyze the spatial and temporal variation of GPP of alpine grasslands on the TP from 1982 to 2015.

Study Area
The TP, located in the southwest of China, is the largest and highest plateau in the world. It covers an area of approximately 2.5 × 10 6 km 2 and has an average elevation of more than 4000 m [30]. At the same time, it is also known as "the Third Pole of the Earth" and "the Roof of the World". Due to its high altitude and complex geographic environment, the TP has formed a unique climate condition characterized by low air temperature, long sunshine duration, and strong solar radiation [31]. Such a special climate environment also promotes the formation of a typical alpine vegetation ecosystem. The dominant vegetation type on the TP is alpine grassland, with an area of about 1.0 × 10 6 km 2 , including alpine steppe and alpine meadow ( Figure 1) [32]. At the same time, the TP is rich in frozen soil, snow, and glacial resources, and their melting water makes the TP the source of several major rivers such as the Yellow River and the Yangtze River. On the other hand, the TP is also extremely sensitive to global climate change. It has been reported that the mean annual temperature on the TP has increased by 0.3 °C per decade over the past 50 years, much higher than other regions of China [5]. Such a sharp warming trend has had a great impact on the alpine vegetation ecosystem, such as advanced spring phenology and enhanced plant carbon fixation [6,23,25].

Meteorological Forcing and Land Surface Initial Conditions for the Biome-BGC Model
In order to estimate the GPP of alpine grasslands at the regional scale, the Biome-BGC model needs spatially explicit forcing data and a vegetation distribution map. The forcing data used to run the Biome-BGC model includes daily meteorological data, land surface initialization data, and ecophysiological parameters. Among them, ecophysiological parameters need to be carefully calibrated to adapt the model to particular applications. For alpine meadows and alpine steppes, both of them are classified as C3 grass in the Biome-BGC model and contains a total of 28 ecophysiological parameters (Table 1). In this study, we used the EC-derived GPP to calibrate these ecophysiological parameters and evaluate the simulation performance of the Biome-BGC model. A detailed description of these data is as follows. In addition, all of the following spatially explicit data was resampled to 5 km × 5 km resolution.

Meteorological Forcing and Land Surface Initial Conditions for the Biome-BGC Model
In order to estimate the GPP of alpine grasslands at the regional scale, the Biome-BGC model needs spatially explicit forcing data and a vegetation distribution map. The forcing data used to run the Biome-BGC model includes daily meteorological data, land surface initialization data, and ecophysiological parameters. Among them, ecophysiological parameters need to be carefully calibrated to adapt the model to particular applications. For alpine meadows and alpine steppes, both of them are classified as C3 grass in the Biome-BGC model and contains a total of 28 ecophysiological parameters (Table 1). In this study, we used the EC-derived GPP to calibrate these ecophysiological parameters and evaluate the simulation performance of the Biome-BGC model. A detailed description of these data is as follows. In addition, all of the following spatially explicit data was resampled to 5 km × 5 km resolution. (1) Daily meteorological data used to drive the Biome-BGC model mainly includes precipitation, temperature, shortwave radiation, day length, and vapor pressure deficit (VPD). They were calculated from the China Meteorological Forcing Dataset (http://westdc.westgis.ac.cn) using the MT-CLIM tool. Among them, the China Meteorological Forcing Dataset has a 0.1 • spatial resolution and a daily temporal resolution, spanning the period from 1982 to 2015, and the MT-CLIM tool is a program released with the Biome-BGC that can be used to estimate missing daily meteorological data that is not available from in situ measurements [33]. In this study, we used the MT-CLIM tool to estimate the day length and VPD that were not included in the China Meteorological Forcing Dataset.
(2) Land surface initialization data mainly contains soil texture (percentage of sand, clay, and silt), soil depth, elevation, shortwave albedo, atmospheric CO 2 concentration, and atmospheric nitrogen deposition. Among them, soil texture and soil depth were obtained from Soil Map Based Harmonized World Soil Database (v1.2) provided by Cold and Arid Regions Sciences Data Center with a spatial resolution of 30 arc-second. Elevation data was acquired from STRM 90m Digital Elevation Database provided by CGIAR Consortium for Spatial Information (available at http://www.cgiar-csi.org). Shortwave albedo data with a spatial resolution of 0.05 • × 0.05 • was obtained from the National Earth System Science Data Sharing Infrastructure, National Science & Technology Infrastructure of China (available at http://www.geodata.cn). Atmospheric CO 2 concentration data was obtained from the National Oceanic and Atmospheric Administration (NOAA) (available at ftp://ftp.cmdl.noaa.gov/ccg/ co2/trends/). Atmospheric nitrogen deposition data at a 0.1 • × 0.1 • resolution was supplied by the Chinese Academy of Agricultural Sciences.
(3) The EC flux data at two flux sites was obtained from the Chinese Flux Observation and Research Network (ChinaFLUX), which was used to calibrate the ecophysiological parameters corresponding to alpine meadows and alpine steppes, respectively. A detailed description of these two flux sites are given in Table 2, and information on the EC measurement system can be found in [34]. The initial measured fluxes (Net Ecosystem Exchange, NEE) were first processed by the ChinaFLUX CO 2 data processing system [35]. Then, the linear interpolation model and the Michaelis-Menten equation were used to fill small gaps (less than 2 h) and large gaps in the time series, respectively [36]. Subsequently, daytime ecosystem respiration was calculated using the Lloyd and Taylor method combined with nighttime observations [37]. Finally, GPP was estimated as the difference between daytime NEE and ecosystem respiration.
(4) The distribution map of alpine grassland was derived from the digitized 1:1 million vegetation map of the TP compiled by Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences ( Figure 1).  Normalized difference vegetation index (NDVI), defined as the difference between the near-infrared and red reflections divided by the sum of the two, is a dimensionless index that can be used to estimate the density of green on an area of land [40]. The Global Inventory Modeling and Mapping Studies (GIMMS) NDVI3g.v1 dataset, with a spatial resolution of 0.083 • × 0.083 • and a temporal resolution of 15-day intervals, has been the most widely used remote sensing product for monitoring vegetation phenology over large areas [41]. It was generated from the Advanced Very High Resolution Radiometer (AVHRR) sensors onboard on several NOAA satellites. At the same time, it was also the latest version of the GIMMS NDVI dataset, spanning the period from January 1982 to December 2015. This dataset has improved data quality at high latitudes, especially in regions where the growing season is shorter than 2 months [42].
Before using the GIMMS NDVI dataset, some preprocessing procedures were needed to reduce the effects of clouds and the atmosphere. Specifically, the maximum value composition (MVC) was first applied to the NDVI dataset to eliminate noises caused by cloud and snow contamination [43]. Then, the Savitzky-Golay filter was used to remove outliers and smooth the NDVI time series [44]. Finally, regions with a multiyear average NDVI less than 0.1 were excluded from the analysis to further minimize the impacts of barren land and non-vegetation [24].

Methodology
The methodology of this study is composed of four key steps. Firstly, the original phenology module of the Biome-BGC was modified by remotely sensed phenology, which was achieved by using the phenological indicators derived from remote sensing to prescribe the phenological states within the Biome-BGC model. Secondly, based on the Biome-BGC model after a phenology module modification, the Morris qualitative sensitivity analysis method was used to screen out those parameters that have negligible effects on the GPP simulation, and then the Sobol' quantitative sensitivity analysis method was used to evaluate the influence of each of the remaining parameters on the simulated GPP and further determine the parameters to be calibrated. Moreover, considering the differences in the physiological characteristics of alpine meadows and alpine steppes, we performed a sensitivity analysis on them separately. Thirdly, we used the simulated annealing optimization algorithm combined with EC-derived GPP data to calibrate those influential parameters selected in the previous step, and set those non-influential parameters to the default values. This calibration process was also performed on alpine meadows and alpine steppes, respectively. In addition, since only two years of EC-derived GPP data was obtained for both alpine meadows and alpine steppes, we used the first year of GPP data to calibrate the influential parameters while using the second year of GPP data to validate the performance of the Biome-BGC model. Finally, we extrapolated the Biome-BGC model after a phenology module modification and parameter calibration to the entire alpine grasslands on the TP to simulate its GPP from 1982 to 2015. Figure 2 shows the organizational flow chart of the methodology of this study, and details will be given in the subsequent sections.

The Biome-BGC Model
Biome-BGC is a process-based biogeochemical model used to simulate the carbon, nitrogen, and water cycles within the terrestrial ecosystems from a single plot to regional and global scales [11]. It uses a daily time-step and was developed from the Forest-BGC model [45]. Previous studies have successfully applied this model to a variety of biomes such as forests and grasslands [17,18,20,46]. There are several critical processes in the Biome-BGC model, including photosynthesis, respiration, evapotranspiration, decomposition, the allocation of photosynthetic assimilates, and mortality [11]. To model the photosynthesis process, the Biome-BGC model converts leaf C into an equivalent leaf area by multiplying the user-defined specific leaf area (SLA) parameters. SLA is a measure of the thickness of a leaf and is defined as the leaf area per unit mass. The Biome-BGC model further partitions leaf C and leaf area into sunlit and shaded fractions using a two-leaf photosynthesis model [47,48] and calculates each photosynthetic process separately. The carbon flux simulated by the Biome-BGC model mainly includes GPP, Net Primary Productivity (NPP), Net Ecosystem Productivity (NEP), and NEE. Among them, GPP is the total carbon fixed by photosynthesis, NPP is the difference between GPP and autotrophic respiration, NEP is the difference between NPP and heterotrophic respiration, and NEE is the difference between NEP and the carbon loss caused by fire. In this study, we used the newly released Biome-BGC v4.2 to simulate the GPP of alpine grasslands on the TP (available at http://www.ntsg.umt.edu/project/biome-bgc.php).

Phenology Module Modification
Vegetation phenology determines the start (SOS) and end (EOS) of the growing season. Accurate representation of phenology is critical for terrestrial biosphere models to simulate the temporal dynamics of biological processes on the land surface [49]. The phenology module used in the Biome-BGC was developed based on the empirical relationship between phenology and meteorological conditions. It defines SOS as the date when both the summed precipitation and the summed soil temperature are greater than the critical values, and defines EOS as the date after July 1st where the average minimum temperature for the next 3 days is less than the annual average minimum temperature [22]. However, this module was found to have significantly delayed the SOS of vegetation in alpine regions, resulting in a severe underestimation of GPP [21]. To this end, we

The Biome-BGC Model
Biome-BGC is a process-based biogeochemical model used to simulate the carbon, nitrogen, and water cycles within the terrestrial ecosystems from a single plot to regional and global scales [11]. It uses a daily time-step and was developed from the Forest-BGC model [45]. Previous studies have successfully applied this model to a variety of biomes such as forests and grasslands [17,18,20,46]. There are several critical processes in the Biome-BGC model, including photosynthesis, respiration, evapotranspiration, decomposition, the allocation of photosynthetic assimilates, and mortality [11]. To model the photosynthesis process, the Biome-BGC model converts leaf C into an equivalent leaf area by multiplying the user-defined specific leaf area (SLA) parameters. SLA is a measure of the thickness of a leaf and is defined as the leaf area per unit mass. The Biome-BGC model further partitions leaf C and leaf area into sunlit and shaded fractions using a two-leaf photosynthesis model [47,48] and calculates each photosynthetic process separately. The carbon flux simulated by the Biome-BGC model mainly includes GPP, Net Primary Productivity (NPP), Net Ecosystem Productivity (NEP), and NEE. Among them, GPP is the total carbon fixed by photosynthesis, NPP is the difference between GPP and autotrophic respiration, NEP is the difference between NPP and heterotrophic respiration, and NEE is the difference between NEP and the carbon loss caused by fire. In this study, we used the newly released Biome-BGC v4.2 to simulate the GPP of alpine grasslands on the TP (available at http://www.ntsg.umt.edu/project/biome-bgc.php).

Phenology Module Modification
Vegetation phenology determines the start (SOS) and end (EOS) of the growing season. Accurate representation of phenology is critical for terrestrial biosphere models to simulate the temporal dynamics of biological processes on the land surface [49]. The phenology module used in the Biome-BGC was developed based on the empirical relationship between phenology and meteorological conditions. It defines SOS as the date when both the summed precipitation and the summed soil temperature are greater than the critical values, and defines EOS as the date after July 1st where the average minimum temperature for the next 3 days is less than the annual average minimum temperature [22]. However, this module was found to have significantly delayed the SOS of vegetation in alpine regions, resulting in a severe underestimation of GPP [21]. To this end, we proposed to incorporate remotely sensed phenology into the Biome-BGC model to better prescribe the phenological states within the model. Specifically, we first extracted the SOS and EOS of alpine grasslands on the TP from 1982 to 2015 based on the GIMMS NDVI dataset, and then we used the SOS and EOS extracted from the NDVI dataset to prescribe the phenological states of vegetation within the Biome-BGC model. Through this modification, the Biome-BGC can effectively overcome the limitations of its original phenology module.
In order to extract phenological parameters from NDVI time series data, several methods have been proposed, such as the NDVI thresholds method [22,50], fitting logistic functions [51], largest NDVI increase method [52], and delayed moving average method [53]. Among them, the dynamic threshold method [50] has been proven to be one of the simplest and most effective methods for the extraction of phenological parameters [54,55]. Therefore, we also used this method to extract the SOS and EOS from the NDVI dataset. The ratio of NDVI is defined as where NDVI ratio represents the output ratio, ranging from 0 to 1, NDVI represents the daily NDVI, NDVI max represents the annual maximum NDVI, and NDVI min represents the annual minimum NDVI in the upward or downward direction, which corresponds to the calculation of SOS or EOS, respectively. Based on the NDVI ratio , the SOS and EOS are defined as the dates when the NDVI ratio reaches the dynamic threshold of 0.2 in the upward and downward directions, respectively [55]. Accordingly, the length of growing season (LOS) can be defined as the difference between EOS and SOS.

Model Parameterization
Despite improvements in the phenology module of the Biome-BGC, additional model parameter calibration is still required to achieve accurate simulation of the GPP. However, the Biome-BGC model typically includes a large number of parameters, and these parameters are generally interdependent, making the calibration process very difficult. Therefore, prior to calibrating the model parameters, it is necessary to identify those parameters that have a significant impact on the simulated variables for further calibration. Global sensitivity analysis can quantify the contribution of uncertainty in model inputs to uncertainty in model outputs and allow for the evaluation of interactions between model inputs [56]. In this study, we also used the global sensitivity analysis to identify those ecophysiological parameters that have a significant impact on the GPP simulation. Then, we used the simulated annealing optimization algorithm in conjunction with EC-derived GPP data to calibrate those selected influential parameters to achieve model calibration.
In the Biome-BGC model, both alpine meadows and alpine steppes are classified as the type of C3 grass, which includes a total of 28 tunable ecophysiological parameters (Table 1). However, it should be noted that there are still some differences in the physiological characteristics between alpine meadows and alpine steppes. Therefore, in order to more accurately quantify the GPP of alpine grasslands on the TP, we performed a parameter sensitivity analysis and calibration on them separately. On the other hand, since only two years of EC-derived GPP data was obtained for each flux station, we used the first year of daily measured GPP to calibrate the influential parameters while using the second year of GPP data to validate its performance.

Sensitivity Analysis
Global sensitivity analysis typically requires a large amount of computational resources and is quite time-consuming. To overcome this limitation, we used a combination of qualitative sensitivity analysis and quantitative sensitivity analysis to identify those parameters that have a significant impact on the GPP simulation. Specifically, we first collected the probability distribution function (PDF) of the ecophysiological parameters of C3 grass from relevant literature, as shown in Table 1 we used the Morris qualitative analysis method [57] to screen out those ecophysiological parameters that have negligible effects on the GPP simulation and set them to the default values. Finally, we used the Sobol' quantitative analysis method [58] to evaluate the influence of each of the remaining parameters on the GPP simulation and further determine which parameters to be calibrated. A detailed description of these steps is as follows.
• Distribution of ecophysiological parameters In addition to the two phenological parameters for prescribing the phenological states within the model, the Biome-BGC includes a total of 28 tunable ecophysiological parameters for C3 grass. In addition, there is no evidence of periodic fires in alpine grasslands of the TP in the past few decades, so we set the fire mortality parameter to zero and did not calibrate it. Finally, a total of 27 parameters are left for the next sensitivity analysis, and their PDFs were collected from relevant literature, as shown in Table 1.
• Qualitative sensitivity analysis based on the Morris method The Morris method is a typical One-at-a-Time sensitivity analysis method [57]. Due to its low computational cost, it is often used as the first step in global sensitivity analysis to remove parameters that have negligible effects on the simulated results. Two sensitivity indices, mean (µ) and standard deviation (σ), are calculated by the Morris method to qualitatively evaluate the relative importance of each parameter. To perform the analysis, the Morris method needs to be executed r × (n + 1) times, where r represents the number of trajectories and n represents the number of input parameters. In this study, r was set to 10 [56] and the value of n was 27. In addition, we used the average of the annual GPP between 1982 and 2015 as the output variable of the Morris method and set a low threshold value of µ = 0.05 to distinguish between influential and non-influential parameters. That is, parameters with µ ≤ 0.05 were considered to have negligible effects on the GPP simulation and were excluded from the next quantitative analysis. Further information on the Morris method can be found in the Supplementary Materials.
• Quantitative sensitivity analysis based on the Sobol' method The Sobol' method is a widely used quantitative global sensitivity analysis method, which can decompose the variance of the model output into the contribution of each input parameter and their interactions [58]. Given that V Y represents the total variance of the model output, it can be decomposed into the following components: where V i represents the variance explained by the ith input parameter, V ij represents the variance explained by the interactions between the ith and jth input parameters, and n represents the number of input parameters. Based on Equation (2), the first order sensitivity index can be defined as S i = V i /V Y , and the higher order sensitivity indices can be defined as ..,n = V 1,2,...,i,...,n /V Y . Accordingly, the total order sensitivity index of the ith parameter can be defined as the sum of its first order sensitivity index and all higher order sensitivity index involving it. In addition, a large difference between the first order sensitivity index and the total order sensitivity index indicates that the parameter mainly affects the model output through interactions. The Sobol' method uses the Monte Carlo [59] sampling scheme to generate random parameter samples. To calculate sensitivity metrics, it requires a parameter set with a sample size of M × (2n + 2), where M represents the number of base samples and n represents the number of input parameters. In this study, M was set to 512. In addition, we used the total order sensitivity index to quantitatively evaluate the effect of each ecophysiological parameter on the GPP simulation and set its threshold to 0.05. That is, parameters with a total order sensitivity index greater than 0.05 were considered as influential parameters and required further calibration.

Parameter Calibration
The simulated annealing algorithm is a heuristic optimization algorithm based on Monte Carlo's iterative solution [60]. It is able to find the optimal parameter values in the parameter space that enable the model to generate the best agreement between the simulation and the observation. In this study, we used this algorithm to calibrate the influential ecophysiological parameters selected by the Sobol' method. Specifically, we first designed a cost function for the simulated annealing algorithm based on the difference between the observed GPP and the simulated GPP, with the following form: where E represents the sum of absolute errors, GPP MODi represents the daily GPP simulated by the Biome-BGC model, GPP OBSi represents the daily GPP observed by the EC technique, and K represents the number of observed values. Subsequently, we imported the PDFs of those influential ecophysiological parameters into the simulated annealing algorithm to generate new parameter values, and further calculated its corresponding cost function value. Then, we judged the cost function value according to the Metropolis criterion to decide whether to accept the new parameter values. Finally, the above two steps were repeated until the termination condition was reached. Through these iterative optimization processes, we have a high probability of obtaining global optimal parameter values.

Biome-BGC Model Simulation
The simulation of the Biome-BGC model consists of two phases, namely the spin-up simulation and the normal simulation. For the spin-up simulation phase, the Biome-BGC model started with constant preindustrial atmospheric CO 2 concentration (286 ppm) and nitrogen deposition (0.0002 kgN/m 2 /yr) and looped through the meteorological data from 1982 to 2015 many times until a steady state was reached [61,62]. For the normal simulation phase, the results of the spin-up simulation were used as initial values for the carbon and nitrogen pools and then normal simulations were performed.
In this study, the original Biome-BGC model without any modifications was called Version-0 (V0), the Biome-BGC model whose phenology module was modified by remotely sensed phenology was called Version-1 (V1), and the Biome-BGC model after a phenology module modification and parameter calibration was called Version-2 (V2). In order to estimate the GPP of alpine grasslands at the regional scale, we first used the EC-derived GPP data to validate the performance of V2. Then, we extrapolated V2 to the entire alpine grasslands on the TP to simulate its GPP from 1982 to 2015.

Remotely Sensed Phenology
The spatial patterns of mean annual SOS, EOS, and LOS extracted from the NDVI time series data were provided in Figure 3. As depicted, the SOS showed a decreasing trend from west to east, consistent with the hydrothermal gradient of the TP. The EOS showed a large spatial heterogeneity across the TP and roughly exhibited a decreasing trend from south to north. Meanwhile, the LOS generally showed a decreasing trend from southeast to northwest, which is also consistent with the hydrothermal gradient on the TP. In addition, according to the statistical results, the SOS was mainly concentrated between the 100th and 150th days of the year, accounting for 95.4% of the alpine grassland area. The EOS mainly occurred between the 290th and 320th days of the year, accounting for 94.7% of the alpine grassland area. Furthermore, the LOS was mainly concentrated between 170 and 210 days, accounting for 85.3% of the alpine grassland area.

Parameter Sensitivity Analysis and Calibration
Based on the Biome-BGC model after a phenology module modification (V1), we used the Morris method and the Sobol' method to identify those parameters that have a significant impact on the GPP simulation for further calibration. The results of Morris sensitivity analysis at Haibei Station and Damxung Station are shown in Figure S1. At Haibei Station, three ecophysiological parameters (Lcel, VPDi, and VPDf) had a Morris index ( ) below the preset threshold of 0.05, and six ecophysiological parameters (Llig, FRcel, SLAsha:sun, Gbl, VPDi, and VPDf) at Damxung Station had a Morris index below the preset threshold. These parameters were considered to have negligible effects on the GPP simulation and were excluded from the next Sobol' sensitivity analysis.
The results of Sobol' sensitivity analysis of the remaining parameters at Haibei Station and Damxung Station are provided in Figure S2. As can be seen, at Haibei Station, the simulation of GPP was most sensitive to C:Nfr, and its large first order and total order sensitivity index values indicate that this parameter affects the GPP simulation either directly or through interaction. The second most influential parameters were C:Nleaf and C:Nlit. Among them, C:Nleaf mainly affects GPP through interaction, while C:Nlit directly controls the simulation of GPP. Interestingly, similar sensitivity analysis results also appeared in Damxung Station, which suggests that there is no significant difference between the parameters that have an important impact on the GPP simulation of alpine meadows and alpine steppes. To further determine the influential parameters to be calibrated, we set the threshold of the total order sensitivity index to 0.05. Through screening, Haibei Station retained seven ecophysiological parameters and Damxung Station retained eight ecophysiological parameters, as shown in Table 3.
Based on these selected influential parameters, we used the simulated annealing optimization algorithm combined with EC-derived GPP data to calibrate these parameters to achieve model parameterization. The calibrated results for Haibei Station and Damxung Station are provided in Table 3.

Parameter Sensitivity Analysis and Calibration
Based on the Biome-BGC model after a phenology module modification (V1), we used the Morris method and the Sobol' method to identify those parameters that have a significant impact on the GPP simulation for further calibration. The results of Morris sensitivity analysis at Haibei Station and Damxung Station are shown in Figure S1. At Haibei Station, three ecophysiological parameters (L cel , VPD i , and VPD f ) had a Morris index (µ) below the preset threshold of 0.05, and six ecophysiological parameters (L lig , FR cel , SLA sha:sun , G bl , VPD i , and VPD f ) at Damxung Station had a Morris index below the preset threshold. These parameters were considered to have negligible effects on the GPP simulation and were excluded from the next Sobol' sensitivity analysis.
The results of Sobol' sensitivity analysis of the remaining parameters at Haibei Station and Damxung Station are provided in Figure S2. As can be seen, at Haibei Station, the simulation of GPP was most sensitive to C:N fr , and its large first order and total order sensitivity index values indicate that this parameter affects the GPP simulation either directly or through interaction. The second most influential parameters were C:N leaf and C:N lit . Among them, C:N leaf mainly affects GPP through interaction, while C:N lit directly controls the simulation of GPP. Interestingly, similar sensitivity analysis results also appeared in Damxung Station, which suggests that there is no significant difference between the parameters that have an important impact on the GPP simulation of alpine meadows and alpine steppes. To further determine the influential parameters to be calibrated, we set the threshold of the total order sensitivity index to 0.05. Through screening, Haibei Station retained seven ecophysiological parameters and Damxung Station retained eight ecophysiological parameters, as shown in Table 3.
Based on these selected influential parameters, we used the simulated annealing optimization algorithm combined with EC-derived GPP data to calibrate these parameters to achieve model parameterization. The calibrated results for Haibei Station and Damxung Station are provided in Table 3.

Site-Level Evaluation
The effects of phenology on GPP simulation for the calibration year and the validation year are shown in Figure 4 and Figure S3, respectively. Among them, the calibration year represents the first year of the acquired EC-derived GPP data, which is used for model calibration, and the validation year represents the second year of the acquired EC-derived GPP data, which is used for model validation. As can be seen from Figure 4 and Figure S3, V0 severely underestimated the GPP values at both Haibei Station and Damxung Station as compared to the EC-derived GPP (marked with red dots). However, the GPP simulated by V1 effectively overcame this limitation, which indicated that remotely sensed phenology can more realistically reflect the phenological characteristics of alpine grasslands on the TP. The reason for the substantial underestimation of GPP in V0 is mainly attributable to the severe lag in the simulation of SOS, which directly leads to a shortened LOS and further affects the carbon fixation process of vegetation. Nevertheless, it should be noted that there was still a large difference between the GPP simulated by V1 and the GPP derived from EC measurements, so further calibration of V1 is still necessary.
year of the acquired EC-derived GPP data, which is used for model calibration, and the validation year represents the second year of the acquired EC-derived GPP data, which is used for model validation. As can be seen from Figure 4 and Figure S3, V0 severely underestimated the GPP values at both Haibei Station and Damxung Station as compared to the EC-derived GPP (marked with red dots). However, the GPP simulated by V1 effectively overcame this limitation, which indicated that remotely sensed phenology can more realistically reflect the phenological characteristics of alpine grasslands on the TP. The reason for the substantial underestimation of GPP in V0 is mainly attributable to the severe lag in the simulation of SOS, which directly leads to a shortened LOS and further affects the carbon fixation process of vegetation. Nevertheless, it should be noted that there was still a large difference between the GPP simulated by V1 and the GPP derived from EC measurements, so further calibration of V1 is still necessary.
A comparison of GPP simulations before and after calibration for the calibration year and the validation year were shown in Figure 4 and Figure S3, respectively. As depicted in Figure 4, V2 significantly improved the simulation accuracy of GPP and showed the best performance in all versions of the Biome-BGC model. Compared with the initial version V0, the RMSE of the final version V2 was reduced by 1.49 gC/m 2 /day at Haibei Station and 0.27 gC/m 2 /day at Damxung Station, and its coefficient of determination 2 ( < 0.01) increased by 0.08 and 0.03 at Haibei Station and Damxung Station, respectively. However, it should be mentioned that the simulation accuracy of V2 is slightly lower than that of V0 at Damxung Station for the validation year ( Figure S3), but it is still much higher than the simulation accuracy of V1. The reason may be attributable primarily to the relatively low carbon fixation rate at Damxung Station, while the GPP simulated by V0 is also underestimated due to the severe lag of the simulated SOS, thus making the GPP simulated by V0 close to the EC-derived GPP.  On the other hand, we also used EC-measured NEE and total ecosystem respiration (TER) to validate the performance of V2. As can be seen from Figure 5 and Figure S4, compared with V0 and V1, V2 significantly improved the simulation accuracy of NEE during the growing season at both Haibei Station and Damxung Station. However, it should be noted that there were still some discrepancies between the NEE simulated by V2 and the NEE measured by EC technique during the non-growing season, which may be attributed to the limitations in the soil respiration module of the Biome-BGC. In addition, as shown in Figure 6 and Figure S5, V2 also greatly improved the simulation   A comparison of GPP simulations before and after calibration for the calibration year and the validation year were shown in Figure 4 and Figure S3, respectively. As depicted in Figure 4, V2 significantly improved the simulation accuracy of GPP and showed the best performance in all versions of the Biome-BGC model. Compared with the initial version V0, the RMSE of the final version V2 was reduced by 1.49 gC/m 2 /day at Haibei Station and 0.27 gC/m 2 /day at Damxung Station, and its coefficient of determination R 2 (P < 0.01) increased by 0.08 and 0.03 at Haibei Station and Damxung Station, respectively. However, it should be mentioned that the simulation accuracy of V2 is slightly lower than that of V0 at Damxung Station for the validation year ( Figure S3), but it is still much higher than the simulation accuracy of V1. The reason may be attributable primarily to the relatively low carbon fixation rate at Damxung Station, while the GPP simulated by V0 is also underestimated due to the severe lag of the simulated SOS, thus making the GPP simulated by V0 close to the EC-derived GPP.
On the other hand, we also used EC-measured NEE and total ecosystem respiration (TER) to validate the performance of V2. As can be seen from Figure 5 and Figure S4, compared with V0 and V1, V2 significantly improved the simulation accuracy of NEE during the growing season at both Haibei Station and Damxung Station. However, it should be noted that there were still some discrepancies between the NEE simulated by V2 and the NEE measured by EC technique during the non-growing season, which may be attributed to the limitations in the soil respiration module of the Biome-BGC. In addition, as shown in Figure 6 and Figure S5, V2 also greatly improved the simulation accuracy of TER, which generated a good agreement with the observed TER. Overall, the Biome-BGC model, after a phenology module modification and parameter calibration, effectively improved the simulation accuracy of carbon fluxes and had a good consistency with the observations, which can be used for GPP simulation of alpine grasslands on the TP.  On the other hand, we also used EC-measured NEE and total ecosystem respiration (TER) to validate the performance of V2. As can be seen from Figure 5 and Figure S4, compared with V0 and V1, V2 significantly improved the simulation accuracy of NEE during the growing season at both Haibei Station and Damxung Station. However, it should be noted that there were still some discrepancies between the NEE simulated by V2 and the NEE measured by EC technique during the non-growing season, which may be attributed to the limitations in the soil respiration module of the Biome-BGC. In addition, as shown in Figure 6 and Figure S5, V2 also greatly improved the simulation accuracy of TER, which generated a good agreement with the observed TER. Overall, the Biome-BGC model, after a phenology module modification and parameter calibration, effectively improved the simulation accuracy of carbon fluxes and had a good consistency with the observations, which can be used for GPP simulation of alpine grasslands on the TP.

Spatial Patterns and Trends of GPP
Based on V2, we simulated the GPP of alpine grasslands on the TP from 1982 to 2015. The spatial pattern of mean annual GPP and its standard deviation (sd) were shown in Figure 7. Generally, the mean annual GPP (Figure 7a) gradually decreased from southeast to northwest, consistent with the

Spatial Patterns and Trends of GPP
Based on V2, we simulated the GPP of alpine grasslands on the TP from 1982 to 2015. The spatial pattern of mean annual GPP and its standard deviation (sd) were shown in Figure 7. Generally, the mean annual GPP (Figure 7a) gradually decreased from southeast to northwest, consistent with the hydrothermal gradient of the TP. At the same time, the sd of GPP (Figure 7b) showed a decreasing trend from east to west, indicating that GPP in the eastern region has experienced greater fluctuations than GPP in the western region in the past 34 years. According to the statistical results, the mean annual GPP of alpine grasslands on the TP was 289.8 gC/m 2 /yr (305.8 TgC/yr), of which alpine meadows and alpine steppes accounted for 461.6 gC/m 2 /yr and 87.9 gC/m 2 /yr, respectively. In addition, most regions (65.3% of the alpine grassland area) had low sd of GPP within 60 gC/m 2 /yr over the study period, mainly distributed in the western region of the TP. Meanwhile, only 12.4% of the alpine grassland areas had sd greater than 120 gC/m 2 /yr, mainly located in the eastern region of the TP.

Spatial Patterns and Trends of GPP
Based on V2, we simulated the GPP of alpine grasslands on the TP from 1982 to 2015. The spatial pattern of mean annual GPP and its standard deviation (sd) were shown in Figure 7. Generally, the mean annual GPP (Figure 7a) gradually decreased from southeast to northwest, consistent with the hydrothermal gradient of the TP. At the same time, the sd of GPP (Figure 7b) showed a decreasing trend from east to west, indicating that GPP in the eastern region has experienced greater fluctuations than GPP in the western region in the past 34 years. According to the statistical results, the mean annual GPP of alpine grasslands on the TP was 289.8 gC/m 2 /yr (305.8 TgC/yr), of which alpine meadows and alpine steppes accounted for 461.6 gC/m 2 /yr and 87.9 gC/m 2 /yr, respectively. In addition, most regions (65.3% of the alpine grassland area) had low sd of GPP within 60 gC/m 2 /yr over the study period, mainly distributed in the western region of the TP. Meanwhile, only 12.4% of the alpine grassland areas had sd greater than 120 gC/m 2 /yr, mainly located in the eastern region of the TP. The spatial distribution of the trends of GPP between 1982 and 2015 was provided in Figure 8. Obviously, the GPP trends on the TP showed a large spatial heterogeneity, which may be mainly attributable to the different changes in climatic conditions. Overall, the mean annual GPP in most The spatial distribution of the trends of GPP between 1982 and 2015 was provided in Figure 8. Obviously, the GPP trends on the TP showed a large spatial heterogeneity, which may be mainly attributable to the different changes in climatic conditions. Overall, the mean annual GPP in most regions (96.5%) showed an increasing trend, with 63.5% of the alpine grassland area increasing significantly. This trend was also reflected in its interannual variations ( Figure 9a). As depicted, the mean annual GPP increased significantly from 1982 to 2015, with a rate of 2.91 gC/m 2 /yr, and it varied from 250.9 gC/m 2 /yr in 1982 to 347.5 gC/m 2 /yr in 2015. Meanwhile, its changing trend was closely related to the changes in precipitation and temperature, as shown in Figure 9b,c. In addition, it should be noted that the mean annual GPP in 1998 had a significant increase as compared to the mean annual GPP of the previous year. This anomaly may be attributed to the EI Nino and La Nina events in 1998. Specifically, the average temperature and precipitation in 1998 were obviously higher than those in 1997, and the SOS of most alpine grasslands was also advanced [23], all of which directly contributed to the significant increase of GPP in 1998. from 250.9 gC/m 2 /yr in 1982 to 347.5 gC/m 2 /yr in 2015. Meanwhile, its changing trend was closely related to the changes in precipitation and temperature, as shown in Figure 9b and Figure 9c. In addition, it should be noted that the mean annual GPP in 1998 had a significant increase as compared to the mean annual GPP of the previous year. This anomaly may be attributed to the EI Nino and La Nina events in 1998. Specifically, the average temperature and precipitation in 1998 were obviously higher than those in 1997, and the SOS of most alpine grasslands was also advanced [23], all of which directly contributed to the significant increase of GPP in 1998.

The Role of Remotely Sensed Phenology in the Biome-BGC Model
Vegetation phenology strongly influences energy, carbon, and nitrogen cycling between the atmosphere and the terrestrial ecosystems. Accurate phenology simulation is a prerequisite for the accurate simulation of carbon fluxes. As mentioned earlier, the original meteorological-based phenology module in the Biome-BGC was found to significantly delay the SOS of alpine grasslands, which seriously affected the simulation accuracy of GPP. This limitation has also been discovered by related to the changes in precipitation and temperature, as shown in Figure 9b and Figure 9c. In addition, it should be noted that the mean annual GPP in 1998 had a significant increase as compared to the mean annual GPP of the previous year. This anomaly may be attributed to the EI Nino and La Nina events in 1998. Specifically, the average temperature and precipitation in 1998 were obviously higher than those in 1997, and the SOS of most alpine grasslands was also advanced [23], all of which directly contributed to the significant increase of GPP in 1998.

The Role of Remotely Sensed Phenology in the Biome-BGC Model
Vegetation phenology strongly influences energy, carbon, and nitrogen cycling between the atmosphere and the terrestrial ecosystems. Accurate phenology simulation is a prerequisite for the accurate simulation of carbon fluxes. As mentioned earlier, the original meteorological-based phenology module in the Biome-BGC was found to significantly delay the SOS of alpine grasslands, which seriously affected the simulation accuracy of GPP. This limitation has also been discovered by

The Role of Remotely Sensed Phenology in the Biome-BGC Model
Vegetation phenology strongly influences energy, carbon, and nitrogen cycling between the atmosphere and the terrestrial ecosystems. Accurate phenology simulation is a prerequisite for the accurate simulation of carbon fluxes. As mentioned earlier, the original meteorological-based phenology module in the Biome-BGC was found to significantly delay the SOS of alpine grasslands, which seriously affected the simulation accuracy of GPP. This limitation has also been discovered by [21], and its reason may be attributed to the fact that the moisture threshold and the temperature threshold in the original meteorological-based phenology module are not applicable to alpine regions like the TP. Specifically, the meteorological-based phenology module in the Biome-BGC assumes that leaf onset occurs when both the summed soil temperature and the summed precipitation are greater than the critical values [22]. However, for the TP, the seasonal thawing of frozen soil and snow cover can provide water for the germination and growth of alpine vegetation, which will make the original moisture threshold become inapplicable. Furthermore, the temperature threshold in the original phenology module of the Biome-BGC was derived from the empirical data for well-suited regions, which may also not well represent the vegetation growth on the TP. Therefore, improvements to the phenology module of the Biome-BGC is required.
Due to the temporal and spatial continuity of remote sensing technology, it has been successfully used for the study of vegetation phenology on the TP [25,26]. In this study, we have used the remotely sensed phenology to prescribe the phenological states within the Biome-BGC model. Results show that the remotely sensed phenology can more accurately represent the phenological characteristics of alpine grasslands than the original meteorological-based phenology module in the Biome-BGC. In addition, different from the previous studies that improved the phenology module of the Biome-BGC [19,20], our study integrated the remotely sensed phenology with the Biome-BGC model, which is more suitable for applications over large areas. This strategy could also be extended to other process-based ecosystem models. However, some limitations of this strategy should also be mentioned. Because remote sensing is an earth observation technique, it can only observe vegetation phenology in a hindcasting manner. Therefore, the Biome-BGC model after a phenology module modification in this study is incapable of making future predictions based on future climate scenarios. At the same time, since NDVI is not sensitive to seasonal changes in most evergreens, this strategy is also not applicable to evergreen species. On the other hand, it is worth mentioning that Sun et al. [63] proposed a multi-factor-driven prognostic phenology module for alpine meadows by combining temperature, precipitation, photoperiod, insolation, and snowfall, which could also be integrated into the Biome-BGC model to overcome the limitations in its original phenology module.

Performance of the Calibrated Biome-BGC Model
Model parameter calibration has long been considered as a difficult task, especially for complex process-based ecosystem models such as the Biome-BGC. In this study, we proposed a practical method for model calibration by combining a global sensitivity analysis with the simulated annealing optimization algorithm. The calibration results show that the carbon fluxes simulated by the calibrated Biome-BGC model are in good agreement with the EC measurements. However, it should be noted that the NEE simulated by the calibrated model is still overestimated during the non-growing season, and its reason can be traced back to the overestimation of TER. Furthermore, the overestimation of TER can be largely attributed to the soil respiration simulation, which is the main source of respired carbon during the non-growing season. The problem of soil respiration module in the Biome-BGC has also been reported in previous studies [17,19,64], and its reason was attributed to either the flaw in the respiration algorithm or the oversimplified simulation of soil temperature [19]. Nevertheless, for the TP, the problem of soil respiration module may also be attributed to the widespread distribution of frozen soil in the non-growing season. Specifically, the presence of frozen soil usually has a large impact on soil temperature, which will further affect the activity of soil microbes and ultimately act on soil respiration. Another study conducted on the TP has also demonstrated the important effects of frozen soil on soil respiration during the non-growing season [65].
On the other hand, some problems regarding model parameter calibration by numerical inversion methods must be acknowledged. Since the mathematical methods used to estimate model parameters usually treat the model as a "black box", that is, the calibration procedure is only based on mathematical optimization criteria without considering the physiological processes within the model, so the calibrated results may not be true [20]. Therefore, some studies suggest avoiding the use of mathematical methods for model calibration and instead using field measurements to determine parameter values [66]. However, for complex process-based ecosystem models like the Biome-BGC, it usually contains a large number of tunable parameters (e.g., 28 ecophysiological parameters for the Biome-BGC), which means that it is almost impossible to obtain all (or even some specific) model parameters through field measurements, especially for applications over large areas. Admittedly, we can choose a subset of parameters for observation, but considering the correlation between model parameters and the large difference in the sensitivity of the parameters to the model output [61], it is still a difficult task, especially when some parameters cannot be observed directly. Therefore, in the absence of alternatives, combining sensitivity analysis with numerical inversion methods for model calibration is still a relatively feasible and reliable approach, which is also the method used in this study.

Uncertainties and Limitations
The spatial pattern of GPP of alpine grasslands on the TP exhibited obvious spatial heterogeneity in our study (Figure 7), consistent with other studies revealed that the GPP of Tibetan alpine grasslands showed a decreasing trend from southeast to northwest [14,16]. In addition, the mean annual GPP of Tibetan alpine grasslands from 2003 to 2008 simulated by [14] using a modified VPM model was 312.3 gC/m 2 /yr, which is slightly lower than the GPP (317.1 gC/m 2 /yr) simulated by our study during the same period, and this small discrepancy may be attributable to the differences in the model used.
On the other hand, some uncertainties still remain in our simulation of the carbon fluxes based on the improved Biome-BGC model. Firstly, since we used the EC-derived GPP data to calibrate the Biome-BGC model, which is essentially obtained through flux partitioning rather than direct observation, the GPP data will inevitably have errors and will further bring uncertainty to the model calibration and GPP simulation. Secondly, disturbances in human activities such as grazing and land use changes can also lead to bias in the GPP simulation. Thirdly, the NDVI dataset used to determine the phenology metrics may also bring uncertainty to the GPP simulation. Specifically, some studies reported that newly launched sensors such as MODIS have better performance than AVHRR for determining phenological parameters in alpine regions [24,67,68]. However, MODIS data was not available until 2000, so we chose to use a single data source to extract phenological parameters to avoid new uncertainties caused by the combination of different sensors. In addition, it is worth mentioning that Kern et al. [69] have improved the quality of GIMMS NDVI3g dataset by applying a local correction to the NDVI3g dataset using MODIS NDVI, which could also be considered in our future work. Finally, uncertainties may also come from the internal structure of the Biome-BGC model. Specifically, due to the widespread distribution of frozen soil on the TP, they will have an important impact on the physiological processes of alpine grasslands, which in turn will affect the carbon cycle and energy flow of alpine grassland ecosystems [25,70]. However, the existence of frozen soil and its effects on the physiological processes have not been considered in the Biome-BGC model, so further improvements to the model structure are still needed. On the other hand, it is worth mentioning that some already improved versions of the Biome-BGC model could also be considered in our future work to reconcile the above issues. For example, Hidy et al. [71] obtained an improved Biome-BGC model (referred to as the Biome-BGCMuSo) by adding a multilayer soil module and a management module to the original model, which could help with the issue of frozen soil and human activities.

Conclusions
In view of the limitations of the Biome-BGC model when applied to alpine regions, we have prescribed its phenological states based upon the remotely sensed phenology and calibrated its parameters to improve the modeling of GPP of alpine grasslands on the Tibetan Plateau (TP). Specifically, we first used the remotely sensed phenology to modify the original meteorological-based phenology module in the Biome-BGC to better prescribe the phenological states within the model, and then we combined the global sensitivity analysis method with the simulated annealing algorithm to calibrate the model parameters. Furthermore, we used the improved Biome-BGC model to simulate the GPP of alpine grasslands on the TP from 1982 to 2015. The results show that the improved Biome-BGC model effectively overcame the limitations in the original model and was capable of reproducing the magnitude and the temporal dynamics of GPP of alpine grasslands. At the same time, the simulated results indicate that the GPP of alpine grasslands on the TP exhibited a large spatial heterogeneity, which generally showed a decreasing trend from southeast to northwest and was consistent with the hydrothermal gradient of the TP. In addition, the mean annual GPP of alpine grasslands on the TP from 1982 to 2015 estimated by our study was 289.8 gC/m 2 /yr (305.8 TgC/yr), and GPP in most regions showed a significant increasing trend. However, it should be noted that although GPP exhibits an increasing trend over the past 34 years, the general trend in net carbon uptake is still unclear given the large uncertainty in soil respiration during the non-growing season. Therefore, in the case that more field-measured data such as soil temperature is available, we also need to consider the effects of frozen