Incorporating water availability into autumn phenological model improved China’s terrestrial gross primary productivity (GPP) simulation

Ecosystem models provide an effective approach to quantify the terrestrial carbon cycle, but the lack of accurate phenological information prevents them from better simulations of the physical processes. Compared with spring phenology (i.e. the start of the growing season, SOS), the vegetation phenology in autumn (the end of the growing season, EOS) is not well-simulated and it is challenging to incorporate vegetation phenology into ecosystem models. The simulation of EOS based on temperature and photoperiod was widely accepted, such as Delpierre et al (2009 Agric. For. Meteorol. 149 938–48)’s model (DM), yet its accuracy has not been fully discussed at a regional scale. Here, we developed a regional autumn phenological model (DMS) with inputs of temperature, photoperiod, and water availability for China’s terrestrial ecosystems. The new DMS model significantly improved the representation of EOS in terms of the lower root mean square error (RMSE), higher model efficiency, and a higher percentage of significant correlation with the referenced EOS. We observed widespread delaying trends of EOS with an average rate of 0.1 d yr−1 for vegetated areas over 2001–2018. We further incorporated the improved EOS into the boreal ecosystem productivity simulator (BEPS) and found that the phenology-modified BEPS model had better performances in predicting annual gross primary productivity (GPP) with ∼28% lower RMSE than the original model when testing against GPP measurements from flux tower sites. From 2001 to 2017, the interannual GPP simulated by the modified BEPS model showed an increasing trend with a rate of 6.0 g C m−2 yr−2. In conclusion, our study proves that water availability is of great significance for modeling autumn phenology, and the incorporation of phenological dates into an ecosystem model is helpful for productivity simulation.


Introduction
Terrestrial ecosystem models provide an efficient approach to quantitively describe the complex interactions between carbon, water, and energy cycles (Liu et al 1997. The ecological carbon sequestration by photosynthesis processes is important for carbon cycle, and accurately estimating the gross primary productivity (GPP) is a key function of the ecosystem models (Gonsamo et al 2013, Liu et al 2018a. An increasing number of studies found that ecosystem productivity is largely controlled by vegetation phenology (Piao et al 2008, Richardson et al 2012. An earlier start of the growing season (SOS) induced by warming tends to accelerate leaf area expanding, increase the nitrogen concentration for the period of plant growth, and facilitate carbon uptake in spring (Richardson et al 2010, Xia et al 2014, Piao et al 2015. Warming in autumn could delay the end of the growing season (EOS) and lengthen plant growth and photosynthesis (Ge et al 2015, Gill et al 2015, Liu et al 2016b, Menzel et al 2020, yet the net carbon sequestration may be partly offset by the substantial carbon consumption through ecosystem respiration in warm autumn (Piao et al 2008, Keenan et al 2014. The timing of vegetation phenology could ultimately alter ecosystem carbon budget, and consequently, an appropriate representation of the plant phenology in ecosystem models is of great significance for the accurate estimation of plant productivity (Richardson et al 2012).
It seemed that the accuracy of ecosystem models was suppressed because the phenological metrics used were oversimple and inaccurate (Richardson et al 2012, Keenan et al 2014. For example, setting thresholds for temperature and radiation (Running and Coughlan 1988), using accumulated temperature , or just following the seasonality of satellite LAI (Liu et al 1997) to determine the start and end of plant growth. These methods usually predicted an earlier onset, a later end, and thus an extended growing season, resulting in an overestimation of productivity in most ecosystem models (Richardson et al 2012, Gonsamo et al 2013. Relative to the overall duration, the exact timing of SOS and EOS were expected to be more significant (Richardson et al 2010, Xu et al 2014. Recently, the simulation of GPP has been improved using new phenological algorithms or remotely sensed phenology (Gonsamo et al 2013, Mao et al 2016. Therefore, it is necessary to integrate accurate phenological information into the ecosystem models.
Long-term observations at sites can provide firsthand detailed records of the phenological events, but such studies are confined to site level. The vegetation phenology extracted from satellite data can provide vast amounts of data on a global scale, however, the longest continuous satellite record can only be traced back to the 1980s. Based on physiological mechanisms and environmental forces, phenological models are capable of rebuilding the incomplete historical phenological data and predicting the potential shifts of the phenological dates under future climate change, which makes the combination of phenological model and ecosystem model more valuable and promising to simulate long-term ecosystem productivity at regional and global scales.
Researchers have made deep explorations of the driving forces of spring phenology and accordingly established spring phenological models (Botta et al 2000, Vitasse et al 2011, Cao et al 2018. Relatively, it is difficult to predict the timing of EOS for the lack of knowledge on the underlying mechanisms (Delpierre et al 2009. White et al (1997) predicted the leaf coloring dates based on a photoperiod-modulated temperature threshold for deciduous broadleaf forests. Dufrêne et al (2005) developed a cold-degree day summation model (CDD) and put it into an ecosystem model. Delpierre et al (2009) modified the CDD model by a function of photoperiod (DM), assuming that leaf coloring is an irreversible process activated by photoperiod. Subsequently, the autumn phenological models based on both temperature and photoperiod were proved to outperform those dependent on only temperature or photoperiod (Caffarra et al 2011, Vitasse et al 2011, Jeong and Medvigy 2014, Lang et al 2019. Most of them were focused on deciduous trees, however, other species or vegetation types were scarcely discussed. Recently, studies have provided solid evidence on other controls on autumn phenology, such as the opposite effects of daytime and nighttime temperature as well as the carryover effect of spring phenology, which thereafter facilitates autumn phenology modeling (Keenan andRichardson 2015, Wu et al 2018). Furthermore, low precipitation or drought often aggravates water stress and advances autumn phenology especially under water-limited conditions (Xie et al 2015, Liu et al 2018b, Peng et al 2019. Several pieces of research proved that incorporating water availability could improve the phenological model in autumn , Xie et al 2018, Chai et al 2019. Similar studies were also carried out in some regions of China, such as the Tibetan Plateau (Chai et al 2019, Lang et al 2019, the Inner Mongolia (Ren et al 2019), temperate regions , and even the whole China (Tao et al 2018). But unfortunately, all of them were based on unevenly distributed meteorological stations, largely constraining the regional modeling and predicting of autumn phenology. A previous study certificated the applicability of a regional-scale phenological model based on satellite-derived phenology (Yang et al 2012). Therefore, other grided products or indices representative of water availability may be feasible for building the regional phenological models (Olsson and Jonsson 2015).
Given the fact that autumn moisture limits the end of photosynthesis  and water availability has cumulative effects on autumn phenology (Peng et al 2019), we hypothesis that by including the water conditions as represented by the standardized precipitation evapotranspiration index (SPEI), EOS could be better simulated and this process would be helpful for the estimation of plant GPP at the regional scale. We, therefore, proposed a new autumn phenological model (DMS) and incorporated it into an ecosystem model (the boreal ecosystem productivity simulator, BEPS) to improve the simulation of ecosystem productivity (i.e. GPP). Overall, our main objectives were: (a) to develop a new autumn phenological model considering the water status; (b) to analyze the spatiotemporal patterns of the predicted EOS across China during the period 2001-2018; (c) to improve the GPP simulation of the BEPS model by incorporating accurate phenological dates; (d) to explore the spatiotemporal patterns of annual GPP simulated by the phenology-modified BEPS model across China over 2001-2017.

Datasets
To test the performance of the phenology-modified BEPS model in GPP simulation, we acquired 72 site-year GPP flux records across 17 flux tower sites in China from the EC_ChinaFLUX_V1.0 dataset (http://chinaflux.org/index.aspx), and the longest record spanned from 2003 to 2017 (table 1). Among the 17 sites, one was for cropland, three sites for forests (i.e. ENF, EBF, and MF), and the rest were for grasslands (figure 1). We also employed a total of 17 155 daily records and 720 monthly observations at 14 of 17 flux sites, excepted for the Haibei-2, Haibei-3, and Hongyuan sites because their detailed data were unavailable.
We used the 16 days MOD13C1 (Collection 6) normalized difference vegetation index (NDVI) dataset from the moderate-resolution imaging spectroradiometer (MODIS) to extract phenological dates (Didan 2015). The 0.05 • product is cloud-free and spans from 2001 to 2018.
The SPEI dataset for 2001-2018 was used as the representative of water availability. It was the normal distribution of the difference between precipitation and evapotranspiration, representing the climatic water balance at multiple temporal scales (Vicente-Serrano et al 2010). The SPEI value with a time scale of i months (i.e. i-month SPEI) indicates the relatively climatic water budget during the continuous i months. For a given pixel, a higher SPEI value normally represents moister climate, while a lower value represents drier conditions. We calculated the SPEI indicator of 1 to 12 month time scales during 2000-2018, with a monthly interval and a pixel size of 0.05 • , using the SPEI package in R language (Vicente-Serrano et al 2010). The climatic data used here were provided by Liu et al (2018a), including the maximum and minimum temperature, precipitation, shortwave radiation, wind speed, air pressure, and humidity.
The BEPS model integrated various data types including remote sensing data, meteorological data, and soil data (Liu et al 1997). The remote sensing data include LAI time series and landcover type data, and the meteorological data include the minimum temperature, maximum temperature, precipitation, shortwave radiation, and humidity at a daily interval. We adopted the 8 days LAI dataset (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) from the Global LAnd Surface Satellite (GLASS) products because it had an overall better prediction in the BEPS model than other datasets (Liu et al 2018a, Xie et al 2019. For the landcover data, we used the MODIS Land Cover Climate Modeling Grid (MCD12C1) Version 6 product which has a spatial resolution of 0.05 • (Friedl and Sulla-Menashe 2015). Other input data from 2001 to 2017 were provided by Liu et al (2018a), and all of them were resampled to 0.05 • .
Three global GPP products were adopted to make further comparisons with our GPP simulation. The first one is the GLASS GPP product with a 0.05 • spatial interval (Liang et al 2013), the second one is the MODIS GPP product from the MOD17A2H Version 6 with a fine pixel size of 500 m (Running et al 2015), and the third one is produced with a revised light use efficiency model (EC-LUE GPP) and has a spatial resolution of 0.05 • (Zheng et al 2020). All of them spanned from 2001 to 2017 with an 8 days interval.

Phenological extraction methods
First, we excluded those pixels with an annual average NDVI of less than 0.1 because they often indicate sparsely vegetated lands. Second, we used a Savitzky-Golay filtering algorithm to reduce the atmospheric noises (e.g. clouds) and smooth the annual NDVI curves (Chen et al 2004). Third, the phenological dates (i.e. SOS and EOS) were derived from the first derivation of the NDVI curve fitted by the double logistic method, and they were also determined as the dates when the NDVI ratio reached 0.5 using the dynamic-threshold approach. More detailed information about these two methods could be found in Wu et al (2018). At last, the mean dates (i.e. SOS MOD and EOS MOD ) calculated by the two methods were used in our later modeling and analyses.

Detecting the influence of water changes on EOS
We compared the correlations between the satellitederived EOS and preseason precipitation as well as between the satellite-derived EOS and preseason SPEI during 2001-2018 using Pearson's correlation analysis. The preseason length of precipitation was defined as i months when EOS had the maximum correlation with the precipitation sum during continuous i months before EOS occurred (including the month EOS occurred). The preseason length of SPEI was directly regarded as the i-month time scale of SPEI at the month EOS occurred. The result was presented in figure S1 (available online at stacks.iop.org/ERL/16/094012/mmedia).

Modeling of autumn phenology
The DM model developed by Delpierre et al (2009) predicts the timing of leaf coloring or senescence using the integration of photoperiod and temperature. The process of leaf coloring or senescence starts from DOY d s when the daily temperature drops to a certain threshold (T b ) and the photoperiod is shorter than the critical length (P s ). The state of leaf coloring or senescence at DOY d, S sen (d), is the accumulated daily rate of leaf coloring or senescence (R sen ) from DOY d s to DOY d, and R sen (d) is a function of temperature and photoperiod at DOY d. When S sen (d) reaches its critical value (Y crit ), then leaf coloring or senescence occurs (i.e. EOS DM ). The whole process could be expressed by the following equations: During the parameter optimization, T b ranged from −2 • C to 36 • C with an interval of 1 • C, while P s started from 9 h to 17.5 h with an interval of 0.5 h.
The parameters x and y could be 0, 1, or 2 whereas they could not be 0 at the same time. The former equation in equation (3) suggests a weakening control of shortened day length (P (d) ) on Y crit while the latter one has the opposite meaning. Here, the set of parameters (i.e. T b , P s , x, y, and Y crit ) were optimal for modeling when the root mean square error (RMSE) between the modeled EOS DM and the satellite-derived EOS MOD was the lowest at the pixel level. We further proposed a new model (DMS) by incorporating the SPEI data into the DM model considering the cumulative effects of water availability on autumn phenology (Peng et al 2019). For a given pixel, the critical value Y crit is linearly linked to the SPEI value at the month of DOY d s during its preseason length In equation (5), the preseason length i was determined when i-month SPEI at the month of DOY d s had the highest correlation coefficient with DOY d s , and it ranged between 1 and 12 months. Similarly, the set of parameters (T b , P s , x, y, a, and b) were determined when the RMSE between the modeled EOS DMS and the satellite-derived EOS MOD was the lowest for each pixel.

Modification of the BEPS model
The BEPS model is a widely applied process-based ecosystem model for water and carbon fluxes simulation (Liu et al 1997, Ju et al 2006, and the module of leaf water potential is the conjunction between water and carbon fluxes (Liu et al 1997). Both the carbon and water cycles in BEPS are regulated by stomatal resistance (Jarvis 1976), which is further regulated by environmental factors including photosynthetically active radiation, atmospheric temperature, vapor pressure deficit, and soil water content. The carbon cycle involves three parts: gross photosynthesis, autotrophic respiration of roots, stems, and leaves, and soil respiration (Liu et al 1997, Gonsamo et al 2013. The key equations of GPP simulation and its relation to soil moisture were presented in the supplementary material. The employment of the BEPS model in China's terrestrial ecosystems has been strictly parameterized using parameter values from literatures and eddy covariance (EC) flux sites measurements (Liu et al 2018a). However, the LAI, an important input of this model, has low values in non-growing seasons, which could lead to positive biases of annual GPP, especially for deciduous forests (Richardson et al 2012). Therefore, the BEPS model requires a better representation of the start and end of the GPP simulation. A previous study confirmed that the overestimation of GPP has been lowered by a new phenological equation, f (p,n) , in Canada's ecosystems (Gonsamo et al 2013). To better predict GPP, we modified this function using the actual phenological dates In equation (6), GPP org is the daily GPP simulated by the original algorithm in the BEPS model, and GPP DMS is the result after modification with f (p,n) . In equation (7), d represents the day of year, T soil represents the daily soil temperature, and T air represents the daily air temperature. In our study, SOS was the satellite-derived SOS MOD , and EOS was the EOS DMS modeled by the DMS model. The BEPS model would run at a daily time interval, and the daily values for a whole year (365 days) were accumulated as the annual GPP.

Evaluation of model performances
The parameter optimization in DM and DMS models was determined when the lowest RMSE between the satellite-derived and modeled EOS occurred. Further, we assessed the accuracy of the DM and DMS models based on three indicators: RMSE, Pearson's correlation (r), and model efficiency (ME): In equations (8)-(10), O i indicates the referenced EOS MOD , P i indicates the modeled EOS, and n is the sum of years from 2001 to 2018. Similarly, the performance of the GPP simulation was evaluated using the RMSE and the determination coefficient (R 2 ) between the GPP flux data (GPP flux ) and those predicted by the original (GPP org ) or phenologymodified BEPS model (GPP DMS ). The predicted GPP for validation was the average GPP of the 3 × 3 pixels where a flux site was located in the central pixel, both of them had the identical landcover type, and surrounding pixels with other landcover types were excluded. Lower RMSE, higher r, and higher ME were associated with better performances of models. Also, a higher r which is statistically significant at p < 0.05 was assumed to be a sign of a well-fitted model for autumn phenology.

Performances of the autumn phenological models
The comparisons of the performances between the DM and DMS models were presented in figure 2. Overall, the EOS DMS modeled by the DMS model had higher r (0.50 vs. 0.37), lower RMSE (6.9 vs. 8.9 days), and higher ME (0.23 vs. 0.12) with the referenced EOS MOD than those metrics for the modeled EOS DM , and these two models showed significant differences among the three metrics using the AVONA analysis (figure 2(a), (c) and (e)). As shown in figure 2( figure 2(e)). Classification into vegetation types showed similar results that EOS DMS had lower RMSE, higher r value, and higher ME for all vegetation types (figures 2(b), (d) and (f)). DBF showed the least changes in the three metrics: r value increased from 0.48 to 0.55, RMSE decreased from 4.5 to 4.3 days, and ME increased from 0.19 to 0.27. Compared to other types, the modeled EOS of DBF had the highest r, lowest RMSE, and highest ME with the referenced EOS for both the DM and DMS models, which is followed by GRA and MF. And the EBF showed the smallest r (0.34 and 0.43), largest RMSE (17.7 and 13.3 days), and smallest ME values (0.10 and 0.18) for both models. Figure 3 presents the distribution of the optimal parameters of the DMS model. In figures 3(a) and (b) lower base temperature (T b ) between −2 • C and 14 • C was found at the TP (26.6%), the range of 14 • C-22 • C was distributed in middle China (17.8%), and the higher values of T b above 22 • C was mostly found in the southeast (17.4%). The critical day length (P s ) increased along the latitude, showing that longer P s of >13 h was mainly observed at the areas of northern 30 • N ( figure 3(b)). For the parameters, x and y (figures 3(c) and (d)), the senescence rate controlled only by photoperiod (x = 0) was observed at 85.8% of the vegetated pixels, the only control from temperature (y = 0) occupied 2.5%, and the remaining areas were controlled by both temperature and photoperiod (x > 0 and y > 0). As shown in figure 3(e), about 65.1% of the areas showed a positive coefficient b. The preseason length of SPEI was centralized at the month of d s occurred and its preceding month, accounting for 52.8% of the vegetated pixels ( figure 3(f)). Figure 4 shows the spatiotemporal patterns of the EOS predicted by the DMS model during the period 2001-2018. Overall, the mean EOS mainly occurred in October (DOY 270-300), accounting for 83.4% of the vegetated lands ( figure 4(a)). Earlier EOS (DOY < 275) was primarily found at higher latitudes, higher altitudes (e.g. TP), and the North China Plain (i.e. NCP). Later EOS (DOY > 290) was observed in the central and warmer southern regions. For vegetation types (figure 4(c)), the earlier EOS was observed for CRO (274 ± 8) and GRA (275 ± 7), while the later dates were found for ENF (289 ± 8) and EBF (291 ± 9) areas.

Simulation of autumn phenology by the DMS model
EOS showed widespread delaying trends during 2000-2018 with a mean speed of 0.1 d yr −1 (0.4 d yr −1 for the significant trend, figure 4(d)). The delayed trend in EOS accounted for 68.7% of the vegetated areas, of which 9.2% were significant. Conversely, about 31.3% of the pixels showed an advanced trend in EOS but only 0.9% was significant. With respect to vegetation types, the significantly delayed EOS occupied the largest proportion for ENF (17.1%), while GRA (5.9%) had the least percentage. Moreover, the EBF and ENF showed faster-changing rates of EOS, 0.26 and 0.25 d yr −1 , respectively, while GRA showed the slowest trend (0.05 d yr −1 ).

Performances of the phenology-modified BEPS model
The comparisons of the simulated GPP using the original BEPS model (GPP org ) and the phenologymodified BEPS model (GPP DMS ) were shown in figure 5. For all sites, the GPP simulated by the phenology-modified BEPS model had a relatively higher R 2 (0.94) with the site-level GPP flux , compared with the original version (0.92). The slope of the regression between the annual GPP flux and GPP DMS was closer to 1.0 (0.99). More importantly, the RMSE between GPP flux and GPP DMS (172.8 g C m −2 yr −1 ) was lower than that between GPP flux and GPP org (239.6), leading to a total reduction of 27.9% in RMSE. In terms of the difference between forests and grasslands, the slope of the linear regression declined from 1.61 to 1.32 for forest sites, while the slope remained almost unchanged at grassland sites. Meanwhile, the RMSE largely reduced with the percentages of 40.7% from the original model (438.5) to the modified version (260.2) for the forest sites, while the reduction was 6.6% at the grassland sites.
For the forest sites, the RMSE decreased from 42.94 to 37.57 g C m −2 month −1 at the monthly scale, and it declined from 1.12 to 1.03 g C m −2 d −1 at the daily scale (figure S2). But it showed little change for the grassland sites at both scales (figure S2). For the cropland sites, the slopes of the regression between the simulated and observed GPP were about 0.77 and 0.88 at daily and monthly scales, respectively, and the RMSE slightly increased at both timescales (figure S2). When comparing to other GPP products at the flux sites, the annual GPP DMS values were relatively higher for cropland and lower for forests than GLASS GPP, and it had higher estimations of GPP at both grassland and cropland sites than the MODIS GPP product ( figure S3(a), (b)). We also found that the annual GPP DMS was generally higher for all three types compared to the EC_LUE GPP product ( figure S3(c)). Figure 6 presents the spatiotemporal patterns of the simulated annual GPP using the phenologymodified BEPS model from 2001 to 2017. The mean annual GPP increased from the northwestern regions to southeastern China (figure 6(a)). Lower values (<300 g C m −2 yr −1 ) were largely located in the north and western regions (e.g. Inner Mongolia and TP) (22.0%), while the areas showing higher values (>1500 g C m −2 yr −1 ) were mostly at some parts of the central and southwestern China, the NCP, and southern parts (17.3%). Compared to other vegetation types (figure 6(c)), EBF produced the most annual GPP (1981 ± 261 g C m −2 yr −1 ), followed by CRO (1386 ± 330), WS (1318 ± 318), DBF (1213 ± 177), MF (1205 ± 300), SAV (1114 ± 297), and ENF (749 ± 298). The lowest annual GPP was found for GRA (434 ± 344).

Simulation of annual GPP by the phenology-modified BEPS model
From 2001 to 2017, we observed an overall increasing trend of GPP with the mean rate of 6.0 g C m −2 yr −2 (11.0 g C m −2 yr −2 for the significant trend, ( figure 6(b)). The annual GPP increased for about 82.1% of the vegetated pixels, and 43.0% were significant (p < 0.05). The remaining showed a declining trend, with only 3.8% being significant. Spatially, higher rates (>20 g C m −2 yr −2 ) were primarily observed in central and southern China. In terms of vegetation types, SAV had the greatest increase of GPP (10.1 g C m −2 yr −2 ), while ENF had the slowest growth rate (3.4 g C m −2 yr −2 ).

Performances of the autumn phenological models 4.1.1. Climatic mechanisms of autumn phenological modeling
Phenology modeling is an irreplaceable way to predict the timing and temporal trend of autumn phenology for future climates. However, the inadequate knowledge of the influential mechanisms on EOS makes it difficult to develop accurate models (Delpierre et al 2009, Keenan and Richardson 2015. Here, considering the crucial climatic effects on autumn phenology, particularly water change, we developed the DMS model based on the temperature, photoperiod, and SPEI data. Our results showed that the DMS model had better performances than the previous DM model developed by Delpierre et al (2009).
Numerous studies supported that using changes in photoperiod and temperature could predict the timing of EOS (Borchert et al 2015, Zhang et al 2015, Fu et al 2018. The seasonal growth of trees tends to cease when the day length drops to a critical value, which could be associated with a series of gene expressions such as the CON-STANS, FLOWERING LOCUS T, LIKE-APETALA1, and BRANCHED1 (Petterle et al 2013, Cubas 2020, Maurya et al 2020. Meanwhile, warmer autumn could delay the timing of leaf senescence, because warming can maintain the activity of photosynthetic enzymes, postpone chlorophyll degradation, and prevent plants from autumn frost (Norby et al 2003, Shi et al 2015, Xie et al 2015. Some studies found   that the growth cessation and dormancy in the Rosaceae family (including apple, pear, mountain ash, and stone fruits) could be induced by only low temperature regardless of the photoperiod (Cooke et al 2012, Maurya andBhalerao 2017). Yet underlying molecular mechanisms need further exploration (Petterle et al 2013). Several studies have proved that an autumn phenological model considering both temperature and photoperiod transcended those with only temperature or photoperiod (Vitasse et al 2011. Water deficiency in summer and autumn could result in earlier leaf senescence in water-limited regions (Liu et al 2018b, Peng et al 2019, which makes it a vital environmental cue for autumn phenology modeling. Vegetation growth responds to water stress mainly through two different hydrological strategies, and it differs among different species (Adams et al 2009). First, some species may accelerate stomatal closure and maintain the high water potential of leaves to reduce water loss, which may result in less photosynthesis and earlier leaf coloring ( However, these highly detailed precipitation data are nearly unavailable beyond meteorological sites, which probably makes the water availability less considered in regional phenology modeling. The SPEI indicator could be an appropriate proxy of water availability as it considers both precipitation and evapotranspiration over a long period. Previous analyses have certified the cumulative and lagged impacts of drought on autumn phenology, especially for those regions without sufficient water supplies (Peng et al 2019). Furthermore, we found quite similar patterns of the correlations between EOS and SPEI/precipitation across China (figure S1), and the percentages of the positive significant correlation were higher for DBF, MF, and non-forested types. These results further suggest that SPEI could be a good representative of water availability and it may improve the phenological modeling in autumn.

Performance of the DMS model
Compared with the DM model, our new DMS model had better performances in predicting autumn phenology given its lower RMSE, higher correlation, and higher model efficiency. The percentages largely improved from 28.9% (DM) to 60.6% (DMS) for the significant correlation between the referenced and modeled EOS. Given the fact that most of the previous models were developed for deciduous trees (Delpierre et al 2009, Jeong and Medvigy 2014, our model considered all main vegetation types across China and the accuracy of the DMS model exceeded the DM model for each type. Generally speaking, DBF had the highest r (0.55), lowest RMSE (4.3 days), and highest ME (0.27). We are aware that the pixel-based DMS model could not match well with those species-specific models due to the mixture of various species within each pixel. Nevertheless, a slightly lower RMSE (4.3 days) was observed for deciduous forests (i.e. DBF) relative to the values of 6-9 days in Jeong and Medvigy (2014), 7-10 days in Liu et al (2019), 9-12 days in Tao et al (2018), and 9-13 days in Delpierre et al (2009). Similarly, the RMSE in grasslands (i.e. GRA) was 6.3 days, which is lower than 11.1 days for alpine meadows . On the contrary, we observed the smallest r (0.43), largest RMSE (13.3 days), and smallest ME (0.18) for EBF, which may be associated with the lower precision of the remotely sensed EOS resulted from the unclear seasonality in canopy greenness.
The distribution of the optimal parameters of the DMS model also provides valuable information (figure 3). We observed a decreasing trend in the base temperature (T b ) with the increment of the altitude across the whole China ( figure 3(a)), suggesting that plants at higher altitudes are more tolerant to lower temperatures. Meanwhile, the P s was relatively longer at higher latitudes ( figure 3(b)), suggesting that the earlier leaf senescence in northern areas may be related to the earlier peaks of CO expression at night and earlier downregulation of FT expression compared to southern regions (Petterle et al 2013, Tedla et al 2020. In terms of the parameters x and y (figures 3(c) and (d)), we found that the senescence rate was controlled only by photoperiod (x = 0) at 85.8% of the vegetated pixels, only by temperature (y = 0) at 2.5% of the areas, and the remaining areas were controlled by both temperature and photoperiod (x > 0 and y > 0). This result may help us understand the dominant impacts of temperature and photoperiod on the process of leaf senescence at a regional scale. For SPEI, the positive coefficient b was found at 65.1% of the vegetated areas ( figure 3(e)), indicating that a moister (or drier) climate would increase (or decrease) the requirement of cold accumulation and consequently postpone (or advance) the date of EOS.

Simulated autumn phenology using the DMS model
The general pattern of EOS predicted by the DMS model was consistent with other studies (Yang et al 2015, Liu et al 2016a, with earlier EOS distributed at higher latitudes (e.g. northeast and northwest), higher altitudes (e.g. the TP), and parts of eastern China (e.g. the NCP). The average delay in EOS was 0.1 d yr −1 across China, and the significant trend was 0.4 d yr −1 , which is in good accordance with the range of 0.12-0.48 d yr −1 over China reported by other studies (Ge et al 2015, Yang et al 2015, Liu et al 2016a. We also found faster delays in EOS for evergreen forests (i.e. EBF and ENF) than for the deciduous trees (i.e. DBF) (Liu et al 2016a, Xu et al 2019. However, the EOS trend for grasslands was the slowest, which is consistent with Liu et al (2016a) but opposite to Ge et al (2015) and Xu et al (2019). Further study should be considered for a better comparison among these results.

Uncertainties of the autumn phenological modeling
Most of the previous models were developed using the ground observations as validations (Delpierre et al 2009, Jeong and Medvigy 2014, and the uneven distribution and limited observational sites prevent them from being applied to a wider spatial range. In contrast, the DMS model referencing remote sensing data enables us to predict EOS and further simulate ecosystem productivity at a regional scale. Yet, more attention should be paid to reduce its uncertainties.
First, the referenced satellite-derived EOS is fundamentally important for autumn phenology modeling. Different datasets and preprocessing methods were proved to have substantial impacts on the variation of EOS (Peng et al 2021). Remote sensing products with finer temporal resolution could provide better accuracy especially for autumn phenology (Wu et al 2017), and the snow/cloud contaminations should be strictly processed (Chen et al 2004, Shen et al 2014. It seems that there exist no best extraction methods for all circumstances . Therefore, it is widely accepted to adopt the averaged value extracted by various methods in order to mitigate the errors of extraction methods among different sites or vegetation types (Liu et al 2016b. Second, we still lack a clear understanding of how changes in temperature, photoperiod, and moisture as well as their interactions affect the exact process of leaf senescence, and these unclear mechanisms limited the accuracy of the process-based DMS model. Third, other appropriate driving factors should also be considered following the latest findings. For example, Borchert et al (2015) implied that the daily insolation (the integration of solar intensity over the day length), instead of day length, should be responsible for plant phenology at all latitudes, which may further improve our phenological modeling.

Performance of the phenology-modified BEPS model
To improve the GPP simulation, we modified the GPP simulation in the BEPS model using the SOS MOD and EOS DMS data and validated the result using the fluxsite GPP. As a result, the phenology-modified BEPS model had better performance in GPP simulation than the original version with the higher R 2 , lower RMSE, and the closer-to-one slope of the regression between GPP flux and GPP DMS (figure 5). For all sites, the modified BEPS model not only had a reduction in RMSE (27.9%) but also reduced the overestimation of annual GPP to some extent given the closerto-one slope (0.99). The results in turn emphasize the important impact of plant phenology on ecosystem carbon uptake.
We also found different accuracy of GPP simulation for different landcover types. The RMSE was largely reduced with a percentage of 40.7% at forest sites, but this percentage was less at grassland sites (6.6%). Similar results were also observed at both monthly and daily scales (figure S2), which probably indicates that adopting exact phenological dates may be more useful and significant for GPP estimation at forest types than grasslands. The slope of the linear regression was above 1.30 at forest sites but closed to 1.0 at grassland sites (figure 5), which may be attributed to the GLASS LAI product used in the BEPS model. The overestimation of GPP was found for plant species with higher tower GPP (>1000 g C m −2 yr −1 ) using the GLASS LAI data, while the underestimation was found for those species with lower tower GPP (<500 g C m −2 yr −1 ) (Liu et al 2018b). For the cropland site, we observed mild increasements in RMSE between the GPP DMS and GPP flux at both the monthly and daily timescales, and the GPP values simulated by both BEPS models were lower compared to the flux data ( figure  S2). The limitations of the BEPS model and mismatch between the satellite observations and harvesting dates may be responsible for these results. On one hand, the BEPS model did not consider lateral movements of soil water  and the effect of farming practices, such as field fertilizer and irrigation (Ju et al 2010). Lacking lateral movements and irrigation would directly influence various stages of the water cycle and the calculation of stomatal resistance, bringing a series of uncertainties to the update of the water cycle and productivity estimation. Fertilizer and irrigation could largely improve crop yield especially in spring and early summer; therefore, the BEPS model may underestimate the GPP of cropland without considering these activities. On the other hand, harvesting would lead to an abrupt drop in GPP, LAI, and NDVI, however, the coarse temporal intervals of 8 days LAI and 16 days NDVI may not capture this signal well (Wu et al 2017), resulting in mismatches between the satellite-retrieved EOS and harvesting time as well as between the GPP flux and simulated GPP by the BEPS model. Therefore, a series of adjustments should be made to improve the BEPS model if researchers intend to estimate productivity at cropland.

Variation in the simulated GPP using the phenology-modified BEPS model
The general spatial distribution of the annual GPP simulated by the modified BEPS model is similar to previous modeling studies: decreasing from southeastern to northwestern China , Yao et al 2018, Liu et al 2018a. The highest annual GPP at EBF and the lowest at GRA also accorded well with previous studies (Li et al 2013, Liu et al 2018a, Han et al 2020. The higher GPP in DBF (1213 ± 177) than ENF (749 ± 298) could be attributed to the former type has higher daily photosynthesis than the latter within the growing season and they own different phenological strategies (Richardson et al 2010). Between 2001 and 2017, the GPP increased at a mean rate of 6.0 g C m −2 yr −2 , which is within the range of 2.48-9.68 g C m −2 yr −2 reported by other studies (Ma et al 2018, Yan et al 2019.
An interesting phenomenon is that the annual GPP was high above 1500 g C m −2 yr −1 in the NCP where the date of EOS occurred quite early. The NCP is one of the main crop production regions in China, and the primary crops include maize, soybean, and cotton in summer and wheat in winter. According to relevant studies (Tao et al 2014, the mean maturity date of summer crops was around DOY 261-262 across the whole NCP, which approximates the modeled EOS DMS for the same region. Although climate warming usually advances the sowing date and shortens the growing seasons, the maturity date of maize, soybean, and cotton postponed with different paces in recent years at most regions of the NCP (Tao et al 2014, Wang et al 2017. These delayed maturity dates of maize and soybean could be associated with farming activities (e.g. deferring the sowing date in spring) and the introduction of cultivars requiring higher thermal accumulation to maintain crop production in response to warming (Tao et al 2014. Therefore, the delayed maturity dates and the prolonged growing periods make a positive contribution to the increase in primary production and crop yields.

Uncertainties of the phenology-modified BEPS model
The EOS data we used here primarily indicate the timing when the greenness of vegetation canopy decreased at the fastest rate during autumn, which is similar to the definitions of autumn phenology in other modelling works (Delpierre et al 2009, Vitasse et al 2011. For example, Delpierre et al (2009) determined the referenced leaf coloring date as the date when over 20%-50% of canopy leaves turn into yellow for 90% of the trees, and Vitasse et al (2011) defined leaf senescence as the date on which 50% of the tree leaves turn to yellow or fall. These autumn phenology metrics maybe not well match with the exact end date of the photosynthesis activity. It is known that leaf senescence is accompanied by the degradation of chlorophyll (Fang et al 1998), and photosynthesis would continue until leaf chlorophyll depleted (Keskitalo et al 2005). The photosynthesis in spring could also start in similar mechanisms. Hence, there may be timing gaps between the start/end date of the actual photosynthesis and the satellite-derived SOS/EOS. Plant photosynthesis during the gaps could be fairly low due to the low temperature, short photoperiod, and lower chlorophyll content, leading to a tiny contribution to the annual GPP. Nevertheless, more accurate phenological dates representing the start and end of photosynthesis should be considered in ecosystem models. Furthermore, we modified the BEPS model using the satellite-derived SOS, and the lack of projected SOS data would prohibit the future prediction of GPP. Thus, employing a robust spring phenological model would be of great significance to the improvement of the ecosystem model and GPP simulation under future climate change. In addition, our study also revealed that the BEPS model should be cautiously applied to cropland given the exclusion of field fertilizing, irrigation, and lateral movement of soil water in the BEPS model and mismatch between harvesting dates and satellite observations as we mentioned above.

Conclusion
In the study, we developed an autumn phenological model and embedded it into the BEPS model to improve the simulation of GPP. We developed the DMS model with temperature, photoperiod, and especially water availability to improve the modeling of autumn phenology. Then, we incorporated the phenological dates (SOS MOD and EOS DMS ) into the BEPS model to improve the ability of GPP simulation. As a result, the new DMS model showed better performances for EOS modeling considering three metrics (i.e. RMSE, r, and ME). We found an overall delayed trend in EOS with an average speed of 0.1 d yr −1 , with 9.2% of the EOS being significantly delayed while only 0.9% significantly advanced for China's terrestrial ecosystems over 2001-2018. Such improvement in EOS was demonstrated to be promising for GPP modeling that annual GPP was better simulated with the phenology-modified BEPS model. We further observed an increasing trend of annual GPP with an average rate of 6.0 g C m −2 yr −2 during the study period. These results emphasize that water availability is an important element to improve the temperature-and photoperiod-regulated phenological model in autumn and the incorporation of accurate phenological information into an ecosystem model is workable and helpful for a better simulation of productivity, which could contribute to a great significance for predicting ecosystem productivity under future climate change.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).