Estimating site-specific optimum air temperature and assessing its effect on the photosynthesis of grasslands in mid- to high-latitudes

The effect of air temperature on photosynthesis is important for the terrestrial carbon cycle. The optimum air temperature for photosynthesis is one of the major parameters in data-driven and process-based photosynthesis models that estimate the gross primary production (GPP) of vegetation under a changing climate. To date, most models use the biome-specific optimum air temperature ( T o p t − b ) parameter. To what degree will the site-specific optimum air temperature ( T o p t − s ) affect GPP simulation results remains unclear. In this study, we estimated T o p t − s by using GPP data from 11 grassland eddy flux tower sites (GPPEC) and satellite vegetation indices (NDVI and EVI). We found that Topt-s parameter values estimated from EVI have good consistency with those from GPPEC at individual sites. We also evaluated the effects of site-specific and biome-specific optimum air temperature parameters on grassland photosynthesis. The results showed that the use of T o p t − s in the Vegetation Photosynthesis Model improved to various degrees in both daily and annual GPP estimates in those grassland flux tower sites. Our results highlight the necessity and potential for the use of T o p t − s in terrestrial GPP models, especially in those situations with large temperature variation (heatwave and cold spill events).


Introduction
The relationship between air temperature and photosynthesis or gross primary production (GPP) of vegetation at the local, regional, and global scales has been studied over many decades (Williams et al 2014;Huang et al 2019). Global warming and climatic extremes (e.g. heatwaves and cold spills) have large impacts on vegetation production across space and time (Mu et al 2011;Jiao et al 2019a;Ryu et al 2019). Accurately quantifying the effects of air temperature on the GPP of vegetation at local, regional, and global scales is critical to improving the modeling of GPP and terrestrial carbon cycles.
In addition to process-based biogeochemical models (Sellers et al 1986;McGuire et al 1995;Zhang et al 2012), satellite-based data-driven biogeochemical models have proven to be a great tool for estimating GPP, as satellite-based sensors provide continuous observations across local, regional, and global scales, especially for regions with limited in situ observations. A number of satellite-based terrestrial carbon models have been developed and used to estimate GPP at various spatial scales (Potter et al 1993;Xiao et al 2004a;Zhao et al 2005;Jiang and Ryu 2016;Ryu et al 2019). Most of these satellite-based GPP models are designed on the basis of the production efficiency concept, also known as light use efficiency (LUE) or the radiation use efficiency concept (table S1, available online at stacks.iop.org/ERL/15/034064/ mmedia). These LUE models estimated daily GPP based on photosynthetically active radiation (PAR) absorbed by vegetation (APAR) and LUE (ε) ( e =ǴPP APAR ). In these models, the LUE (ε) parameter is often estimated as the product of the potential or maximum LUE (ε 0 ) and a few down-regulation scalars such as temperature, precipitation, soil moisture, vapor pressure deficit, and leaf age (Monteith 1972(Monteith , 1977Prince and Goward 1995). Temperature constraint, represented by temperature scalar (T s ), has been used in most LUE models (table S1). The typical response of leaf-level photosynthesis to temperature can be described by a bellshaped relationship (Berry and Bjorkman 1980;Cox 2001;Clark et al 2011;Landsberg et al 2011). Because of the limitation of electric transport and Rubisco activity, plants tend to have low photosynthesis at cool temperatures, and increase to a maximum rate at optimal temperatures in the 20°C-30°C range and then decrease again at very high temperatures (Fitter and Hay 2012). This optimum temperature (T opt ) response has been described across a wide range of plant species (Kirschbaum and Farquhar 1984;Battaglia et al 1996;Fitter and Hay 2012), and ecosystem biomes (Huang et al 2019). In most satellite-based LUE GPP models, T s has been defined as a function of T opt , maximum temperature, and minimum temperature for vegetation growth. As reported, the temperature scalar (T s ) is more sensitive and is more highly governed by the choice of T opt than by maximum and minimum temperatures (Raich et al 1991;Zhang et al 2017b). Thus, the choice of T opt largely affects the T s and finally affects the accuracy of GPP estimates in the models.
Biome-specific optimum air temperature parameters (T opt−b ) have been used to calculate biomespecific temperature scalars (T s−b ) in biogeochemical models (table S1) (King et al 2011), and GLO-PEM (Prince and Goward 1995). Several studies have reported that the use of biome-specific parameters introduced an inaccurate derivation of ε and could be one of the potential error sources in GPP data product (Heinsch et al 2006;Sjöström et al 2013). Note that a vegetation biome often covers a very large geographical domain, and vegetation may adapt to its local climate over years and thus develop a sitespecific optimum air temperature (T opt−s ) (Huang et al 2019). There is a need to quantify the range of T opt−s parameter values and the difference between the T opt−s and T opt−b parameters. There is also a need to quantify the potential of using the T opt−s parameter to improve GPP estimates in the models. To date, only a few studies have reported the use of T opt−s in estimating GPP (Potter et al 2003, Sasai et al 2005. However, these studies have not used the GPP estimates and air temperature from the eddy flux tower sites to quantify the relationship between air temperature and GPP and estimate T opt−s parameter values. Therefore, our knowledge of the T opt−s parameter values and the potential of using the T opt−s parameter to improve GPP estimates is still very limited. A number of scientific questions need to be addressed: what is the most appropriate method for estimating T opt−s from both GPP data in the eddy flux tower sites (GPP EC ) and satellite datasets (NDVI, EVI)? What are the differences between site-specific temperature scalar (T s−s ) and biome-specific temperature scalar (T s−b ), and to what degree does T opt−s affect GPP estimates in the data-driven models?The answers to these questions will help improve our understanding of many aspects of terrestrial carbon cycling, such as the impacts of climatic extremes (e.g. heatwave, cold spill) on the seasonal dynamics and inter-annual variation of GPP.
Grasslands in mid-to high-latitude regions are sensitive and vulnerable to climate variability, and temperature is a major climate factor controlling GPP (Yi et al 2010). Also grasslands have the largest inter-annual variation of gross and net primary production among the major ecosystem types (Fridley et al 2016;Hufkens et al 2016;Knapp et al 2017). Grasslands in these regions have high uncertainties in satellite-based GPP estimates. Compared with in situ flux observations, studies have found that the MOD17 GPP algorithm underestimated grassland GPP from sites to regions Zhu et al 2016;Zhu et al 2018). The VPM GPP product added a C3/C4 ratio for the parameter ε calculation and significantly improved grassland GPP estimates (Zhang et al 2017b). However, larger uncertainties still exist in mid-to high-latitude grassland VPM GPP estimates . The large uncertainties in grassland GPP estimates directly hinder our understanding of inner-and inter-annual GPP dynamics, and affect our assessment of ecosystem response to climate variability. For example, an analysis using MOD17 GPP showed large carbon losses for the US in 2012 because of the warm spring and dry summer (Wolf et al 2016), while the VPM GPP showed a slight carbon uptake . In this study, first we quantified T opt−s parameter values in 11 grassland sites in mid-to high-latitude regions, and compared the T opt−s and T opt−b parameters. Our hypothesis is that the T opt−s parameter for photosynthesis of mid-to high-latitude grasslands varies among the sites and differs substantially from the commonly used T opt−b . In order to explore the effects of the methods that are used to estimate T opt−s , and identify potential data sources for T opt−s calculation across the globe, T opt−s values were calculated and compared with multiple data sources (GPP EC , MODIS NDVI, and EVI) and different methods. Second, we assessed the effects of the T opt−s parameter on GPP estimates in these grasslands sites. Our hypothesis is that the T opt−b parameter value may result in a large overestimation or underestimation of GPP of grasslands in previous GPP products, depending upon the differences between T opt−s at individual sites and T opt−b . The VPM, which was developed under the LUE concept and satellite datasets (Xiao et al 2004b;Zhang et al 2017b), was used to estimate GPP for each site over several years. We also estimated and compared VPM GPP products driven by the two different types of T , s including both -T s s (GPP s ) and T s−b (GPP b ). The results from this study may help improve T opt−s parameter estimates and GPP estimates in the grasslands.

Materials and methods
2.1. Study sites Data from grassland flux tower sites at mid-to highlatitudes were used in this study, and the details for these sites are described in the FLUXNET-2015 dataset. We selected the flux sites based on the following criteria: (1) the site has obvious seasonal changes, winter (daily daytime mean temperature (T DT ) lower than 0°C) lasting at least 2 months for each year; (2) land cover type at the site is homogeneous within the MOD09A1 (500 m) pixel (figure S1); (3) the site has had continuous observation for at least 1 year. In this study, we selected and analyzed 11 grassland sites. Spatial distribution and meteorological information of all the flux tower sites used in the analysis are shown in figure S2 and table S2.

Meteorological data and GPP data from the flux tower sites
The FLUXNET-2015 dataset provides meteorological data, water flux, and CO 2 flux data at half-hourly, hourly, daily, and yearly intervals. We visually checked the tower observations, and the values with low quality such as those with the same values in a whole year were removed. We also calculated daily downward surface solar shortwave radiation (  s s ), daily daytime mean temperature (T DT ), and daily GPP (GPP EC ) which were calculated with the variable USTAR filtering approach and daytime portioning method (Kumar et al 2016). Then, 8 day  s , s T DT and GPP EC were generated from daily products respectively, and used in the VPM GPP simulation and comparison.

MODIS vegetation indices
This study used the MODIS land reflectance product MOD09A1 V006 (500 m spatial resolution and 8 day intervals) (Vermote 2015).
To filter out poor quality observations, we firstly identified those affected by ice, snow, and clouds using the quality control (QC) layer ( 2.4. Methods for estimating site-specific optimum temperature for GPP Biome-specific optimum air temperature (T opt−b ) was used at 27°C as reported in the global VPM GPP product (Zhang et al 2017b). Site-specific optimum air temperature (T opt−s ) was estimated from the analyses of temperature, GPP EC , EVI, and NDVI data at individual flux tower sites. We developed two new methods to estimate T opt−s , namely, the 95% maximum method and the generalized additive model (GAM) regression method. In order to make a comparison with a previous study (Potter et al 2003), we also estimated the T opt−s from NDVI following the method in the CASA model, which is denoted as ---T . opt s CASA NDVI In the response curve between the daily air temperature (x-axis) and GPP or vegetation indices (y-axis) (figure 1), we define the site-specific optimum air temperature as the daily air temperature when GPP or vegetation indices reach their peak value within the growing season.
With the 95% maximum method, we firstly found the maximum values of GPP EC (GPP EC-max ), or EVI (EVI max ), for each site. We calculated the optimum temperature as the daily daytime mean temperature (T DT ) during those observations with GPP or EVI values equal to or higher than 95% GPP EC-max or EVI max (figures 1(a) and (b)). Estimated T opt−s using the 95% maximum method from GPP EC and EVI are denoted as T opt s CASA NDVI was defined as the average monthly T DT when GPP EC-max or EVI max occurred ( figure 1(c)).
With the GAM regression method, the relationship between the GPP EC values (or EVI values) and the T DT at a site over all the years were determined using a cyclic penalized cubic regression spline smooth model in R software. The optimum temperature for this site was then defined as the T DT when GPP EC (or EVI) reached the maximum value in the GAM regression line (figures 1(d), (e)). T opt−s estimated by the GAM method from GPP EC and EVI are denoted as . In the VPM model, daily GPP is estimated by APAR by chlorophyll in the canopy (APAR chl ; APAR chl =FPAR chl * PAR) and LUE (ε), see equations (4) and (5): where ε is LUE, FPAR chl is the fraction of PAR absorbed by chlorophyll, and PAR is the photosynthetic active radiation. EVI is used to estimate FPAR .
chl Temperature stress (T s ) and water stress (W s ) are used to downscale maximum LUE (ε 0 ) and estimate ε.
T s is calculated using the temperature response equation documented in the Terrestrial Ecosystem Model (Raich et al (1991)), as shown in equation (6): T  T T  T T  T T  T T  T T  6   s   min  max   min  max  opt  2 where T is the daily daytime mean air temperature (°C); T , min T , max and T opt are the minimum, maximum, and maximum air temperatures for photosynthesis, respectively. The biome-specific parameters used in the global VPM GPP simulations came from the biome-specific look-up table, and the T opt for grasslands was set as 27°C (Zhang et al 2017b). Four groups of site-specific T s ( -T s s ) were calculated using the two methods (95% max and GAM) from GPP EC and EVI, and are denoted as

Estimation of site-specific optimum air temperature from GPP EC and vegetation index
We estimated T opt−s for each site with the three methods using GPP EC , EVI, and NDVI (table S2). The results ( figure 2) showed that the T opt−s values showed a large difference within the grassland sites, and the estimates of T opt−s for individual sites were very different form the T opt−b used in the global VPM GPP product (27°C). For the estimates of T opt−s based on every method, the difference between the highest and lowest T opt−s of the 11 grassland sites was larger than 10°C . T opt−s calculated from EVI and NDVI were significantly correlated with T opt−s from GPP EC when using same estimation method (root mean square error (RMSE) values are from 1.58-3.28°C). As shown by the linear regression results (RMSE, R 2 , P-value), T opt−s estimates from NDVI using the two methods developed in our study,  were more consistent with T opt−s estimates from GPP EC ( --T opt s GPP EC ) than those from the CASA model ( ---T opt s CASA NDVI ) (figures 2(c) and (d)).   T opt s GAM NDVI are the T opt-s from eddy covariance GPP EC , EVI, and NDVI using the GAM regression method.

Discussion
T opt was generally studied and estimated along the level of organization of species, community, and ecosystem. The studies indicated that T opt varies across species and across ecosystems (biomes) (Kattge and Knorr 2007;Lin et al 2012), and T opt−b was used in the biogeochemical models. Different from most previous studies, our study explored and discussed the variability of T opt across sites within a biome. Our results showed large differences of T opt across sites within a biome, and thus supported the urgent need to address T opt−s in a global terrestrial ecosystem study. In addition, the ecosystem-level T opt−b parameters in previous global process-based ecosystem models were directly scaled from the leaf-level T opt−b parameters, in which the T opt−b values at the ecosystem level were found to be consistently lower than those at the leaf level and varied spatially (Huang et al 2019). Our study introduced the methods by using satellite datasets for ecosystem-level T opt−s extraction. The new methods provide a new way and results for future ecosystem T opt studies. Previous studies have suggested gradually changed T opt values along the latitude, while our study did not find a clear relationship between T opt and latitude, annual precipitation, and temperature for the 11 grassland sites (figure S3). This is likely caused by the limited number of grassland sites, or due to grasslands being sensitive to both temperature anomalies and water supply and cannot be well explained by a single climate factor (Hufkens et al 2016;Green et al 2019).
In recent years, many approaches have been developed to reduce the impacts from biome-specific lookup table parameters and coarse image resolutions in GPP estimates, such as readjusting biome-specific parameters (Sjöström et  Our study contributed the LUE estimates by adjusting the temperature parameter and therefore temperature scalars, which was a less considered direction. Even though the CASA model has already tried to use T opt−s instead of T opt−b in the Net Photosynthesis Productivity (NPP) products (Field et al 1995), the two methods (95% max and GAM regression) developed in our study improved the estimates of T opt−s , significantly. Compared with the T opt−s estimated from NDVI in the CASA model, T opt−s estimated from EVI was more consistent to T opt−s from GPP EC (figure 2), which indicates that EVI is a reliable indicator for T opt−s estimation in space, which could contribute to a , and GPP EC . Simple linear regression models were used at each eddy covariance site, and R 2 and RMSE (g C/m 2 /day) were shown. *** means a P-value less than 0.001. large-scale GPP simulation in the future. Because T opt−s in the CASA model has usually been defined as the monthly mean temperature when NDVI reaches its maximum (Yan et al 2015), thus ---T opt s CASA NDVI had more errors than that with a 95% max and GAM regression (figure 1). What is more, NDVI was more affected than EVI especially at regions mixed with complicated background information (Chang et al 2019). As a previous validation study has proved that the global GPP VPM product with T opt−b has been more reliable than GPP CASA when compared with GOME-2 SIF data, our GPP VPM with T opt−s could be much more competitive in model comparison studies . It is important to apply the T opt−s estimation methods in other land cover types, and explore the effects on GPP simulation. Both the datasets and methods in this study have widely applicability in other land cover types.

SiteID
Accurate T opt−s estimation is a reasonably reliable way for improving GPP estimates. A CASA model research study improved NPP by about 50 g C m −2 yr −1 at China's Shennongjia Forestry District in the Hubei province by slightly improving the T opt−s estimation method, in which the T opt−s was defined as the mean temperature during the period of mature stability . Our results indicated that using T opt−b in previous VPM GPP studies could lead to an underestimation of GPP of 25% for grassland ecosystems annually (figure 4(e)). But we found that even though the use of T opt−s improved GPP estimation and resulted in higher GPP values than using T opt−b in most grassland sites, GPP VPM with T opt−s was still lower than GPP EC from eddy covariance observation for many of the 8 day intervals (figures S4(a)-(d)), and -GPP VPM s was about 5%-12% lower than GPP EC annually (figure 4(e)). The annual underestimation mostly occurred in the higher GPP years with 1400 g C m −2 yr −1 at AT-Neu (2002 and CH-Oe1 (2002, which could be caused by the inter-annual and inner-annual variability of C3/C4 composition which are not well recognized in the models Zhu et al 2018). Åt AT-Neu (figure S6) and CH-Oe1 (figure S7), the start of the season and end of the season from GPP EC and GPP VPM agrees well with each other, but the magnitude differs substantially between them within a few years (e.g. 2002, 2003, 2004. Both shortwave radiation data and vegetation index data do not support a very high GPP EC during the 8 day periods of those years. We used daily GPP portioned by net ecosystem exchange (NEE) in the flux tower sites which has been reported to have errors or uncertainties in some observations (Reichstein et al 2005). Here, we would like to attribute the quality of GPP EC data as a major source of the large discrepancy between annual GPP VPM and GPP EC in some years. The daily GPP data showed that the abnormal GPP EC values could be caused by the intensive rainfall (figures S8 and S9). The consistency between GPP EC and climate data and remote sensing data is important for us to evaluate GPP EC data. However, the use of T opt−s was slightly overestimated for the years with lower annual GPP. The overestimation for low GPP years, mostly occurred at IT-Tor (2009-2013, and could be related to the water stress or lower annual precipitation in these years (628-818 mm yr −1 ) relative to the multiyear mean annual precipitation (920 mm yr −1 ). Under drought conditions, there could be a lower T opt−s than in normal years. Further studies are needed to explore the possible ways to improve GPP estimation at the ecosystem scale. Other likely sources of uncertainty in data-driven GPP products include for example the model structure , meteorological input datasets (Anav et al 2015), and seasonal dynamic of LUE (Wei et al 2017). Many novel approaches have been developed to reduce uncertainties in GPP estimates. For example, a study estimated GPP by only using PAR and EVI (Ma et al 2014). The Photochemical Reflectance Index was found to be significantly correlated to LUE, and was effective in detecting seasonal carbon fluxes in evergreen ecosystems where FPAR and greenness-related vegetation indices change little (Garbulsky et al 2011;Middleton et al 2016). NIRv was better correlated to modeled MODIS FPAR than NDVI and significantly correlated to GPP (Badgley et al 2017), and has been used for GPP estimates globally in 0.5° (Badgley et al 2018). Also, significant linear relationships between GPP and OCO-2based SIF product (GOSIF) contributed to the work that estimated GPP in 0.05°using GOSIF (Li and Xiao 2019). Further studies are needed to explore the possible ways to improve GPP estimation at the ecosystem scale.
The satellite-based -GPP VPM s product with higher estimate accuracies could be more reliable for studying the impacts of climate variability, especially extreme climate events, on the ecosystem. Here, we take drought, which is expected to show an intensified frequency and consequences under climate change (Jiao et al 2016;Jiao et al 2019b), as an example for discussing the possible contributions of our study in a future study. Previous studies based on three different global GPP products reported that the impact of drought on terrestrial primary production was underestimated by satellite-based LUE GPP models (Turner et al 2005;Mu et al 2007;Sims et al 2008). The reason for the underestimation is that these GPP models did not simulate the water balance, or did not account for the direct effects of soil moisture in addition to VPD and changes in greenness (Jiao et al 2019a;Stocker et al 2019). Our study found that GPP VPM computed with T opt−s for the years with higher precipitation showed a greater improvement than for the years with lower precipitation ( figure S5). This result indicated that the T opt−b used in previous global GPP simulations might finally underestimate the decrease of GPP from a normal year to a drought year, which could be one of the reasons for the underestimation of drought impacts on ecosystem productivity. As known, when drought occurs, it is often accompanied by higher temperature (Zhang et al 2017a). The plants thus actually suffer both water stress and temperature stress under drought. As the drought condition T opt−s was different with and lower than T opt−b , the use of T opt−b might not capture well the effect of increasing temperature on GPP, and therefore resulted in a greater underestimation. Future GPP models need to consider the comprehensive impacts from multi-parameters such as temperature, water, canopy structural, leaf nitrogen, and chlorophyll content.

Conclusions
Our study explored the estimates of T opt−s using a satellite and the potential of using T opt−s in estimating the GPP of grasslands. We found that EVI has a similar performance with in situ measured GPP EC for determining photosynthesis T opt−s . We also compared the differences in T opt−s values using different extraction methods and different data sources. Our results provide references with data sources and methods for reliable T opt−s estimation and more accurate GPP simulations at the site and global scales. T opt−s values differ among sites and differ from T opt−b significantly. We found a significant improvement in the accuracy of GPP estimates for grasslands by using T opt−s rather than T opt−b . We suggest that terrestrial ecosystem models should account for site-specific temperature parameters. As the climatic impacts on ecosystems have always been assessed by GPP anomalies, an improved GPP product would help us better understand the impacts of extreme events on terrestrial ecosystem carbon cycles, and better manage terrestrial ecosystems.

Data availability
Any data that support the findings of this study are included and discussed within the article.