Non-uniform seasonal warming regulates vegetation greening and atmospheric CO2 amplification over northern lands

The enhanced vegetation growth by climate warming plays a pivotal role in amplifying the seasonal cycle of atmospheric CO2 at northern lands (>50° N) since 1960s. However, the correlation between vegetation growth, temperature and seasonal amplitude of atmospheric CO2 concentration have become elusive with the slowed increasing trend of vegetation growth and weakened temperature control on CO2 uptake since late 1990s. Here, based on in situ atmospheric CO2 concentration records from the Barrow observatory site, we found a slowdown in the increasing trend of the atmospheric CO2 amplitude from 1990s to mid-2000s. This phenomenon was associated with the paused decrease in the minimum CO2 concentration ([CO2]min), which was significantly correlated with the slowdown of vegetation greening and growing-season length extension. We then showed that both the vegetation greenness and growing-season length were positively correlated with spring but not autumn temperature over the northern lands. Furthermore, such asymmetric dependences of vegetation growth upon spring and autumn temperature cannot be captured by the state-of-art terrestrial biosphere models. These findings indicate that the responses of vegetation growth to spring and autumn warming are asymmetric, and highlight the need of improving autumn phenology in the models for predicting seasonal cycle of atmospheric CO2 concentration.


Introduction
Temporal dynamics of atmospheric CO 2 concentration, climate and terrestrial carbon (C) cycle are strongly linked in the present (Schneising et al 2014) and past (Montañez et al 2016) Earth systems. For example, the recent inter-annual variability of atmospheric CO 2 growth rate is largely caused by fluctuations in terrestrial CO 2 uptake (Myneni et al 1997, Keenan et al 2016, which is mainly driven by variations in climate (Poulter et  Recently, non-uniform warming trends among seasons have been detected over the northern lands (Xu et al 2013, Xia et al 2014. Given that climate warming in different seasons would influence vegetation growth differently (Xu et al 2013, Xia et al 2014, Cai et al 2016, the role of seasonal non-uniform warming in affecting the vegetation growth as well as the recent changes of the atmospheric CO 2 amplitude remains unclear. Some recent evidence has implied weakening correlations of MAT with vegetation growth and atmospheric CO 2 concentration in the past three decades. First, the MAT across the northern high latitudes has kept rising whereas vegetation greenness has begun to decline since late 1990s (Bhatt et al 2013, Jeong et al 2013. Second, the advanced spring phenology in response to climate warming has been reported to diminish at northern high latitudes over the last two decades (Fu et al 2015, Wang et al 2015). Third, a weakening inter-annual correlation of temperature with vegetation greenness  or spring ecosystem CO 2 uptake (Piao et al 2017b) has been detected at northern latitudes during recent years. In northern temperate ecosystems, a negative correlation between the atmospheric CO 2 amplitude and temperature anomalies during 2000s has been found through the analysis of space-borne atmospheric CO 2 measurements (Schneising et al 2014). Thus, it is important to examine how the temperature changes in different seasons have contributed to such weakening correlations of MAT with vegetation growth and seasonal CO 2 amplitude in recent years.
In this study, we investigated the relationships between changes in the seasonal CO 2 amplitude, vegetation greenness and seasonal air temperature in northern lands (>50°N) during the last three decades. The analyses were based on long-term monitoring records of atmospheric CO 2 , global gridded climate datasets, and satellite-derived Normalized Difference Vegetation Index (NDVI). We also examined the relationships between temperature and gross primary productivity (GPP) in five terrestrial biosphere models (TBMs), which have been commonly incorporated into Earth system models for future projections of climate and atmospheric changes.
The in situ [CO 2 ] observations at the BRW site were collected hourly. The daily and monthly data provided by NOAA were averaged from the hourly observations. This study used the monthly data to derive the CO 2 amplitude. To avoid the biases of various processing algorithms, we collected the estimates of [CO 2 ] amplitude from Forkel et al (2016) and the GLOBALVIEW products (figure S2). In Forkel et al (2016), the time series of daily [CO 2 ] records were first fitted with polynomial and harmonics functions and then de-trended with the Fast Fourier Transformation method (Thoning et al 1989, Thoning et al 2015. In each calendar year, the [CO 2 ] amplitude was calculated as the peak-to-trough difference of the de-trended seasonal cycle. The data in the GLOBALVIEW-CO 2 product has already been smoothed, interpolated and extrapolated (Masane and Tans 1955). GLOBALVIEW-CO 2 provides observations at 7 d intervals. We obtained its [CO 2 ] amplitude as the difference between the maximum and minimum weekly CO 2 data in each calendar year. According to the Theil-Sen analysis, the long-term trends of [CO 2 ] amplitude were consistent among these processing methods (figure S2).

Temperature analysis
The temperature trends were analyzed based on the latest version of the CRU temperature data (CRU TS4.0). It is gridded with a spatial resolution of 0.5°×0.5°at a monthly time step. This product is gridded using the Angular-distance weighting interpolation (Harris and Jones 2017) based on the observations collected from 2600 stations worldwide (Harris et al 2014). The CRU climate products have been widely used for phenology analysis and for driving different types of ecosystem models (Koven et al 2011, Fu et al 2014, Xia et al 2017. In this study, seasonal temperature was averaged from monthly temperature following the definition of four seasons: Spring, March-May; Summer, June-August; Autumn, September-November; Winter, December-February. Besides, MAT of a given year was defined as the average of the monthly temperature from January to December. The Theil-Sen estimator and Mann-Kendall trend test were applied in detecting the temporal trends of the seasonal temperature (see more details in section 2.6).

Satellite derived NDVI
The normalized difference vegetation index (NDVI) is widely used as an indicator for vegetation productivity (Myneni et al 1997, Zhou et al 2001. It is calculated as the normalized ratio between near infrared and red reflectance bands (Tucker 1979, Tucker et al 2005. The NDVI used in this study is from the Advanced Very High Resolution Radiometer sensors, which has the longest record of continuous satellite data since 1981. Here, we used the newest version of GIMMS NDVI dataset (NDVI3g) (Tucker et al 2005, Pinzon andTucker 2014). It is a global product at spatial resolution of ∼8×8 km and temporal resolution of 15 d. The NDVI3g has been widely used for analyzing vegetation changes in recent years (Tucker et al 2005, Peng et al 2013, Wang et al 2014. The maximum value composites method (Holben 1986) was applied to merge segmented data strips to half-monthly values. To lessen the impacts of sparse soils and snows on vegetation, following Zhang et al (2013), the areas with multiyear average NDVI less than 0.1 in the northern lands (>50°N) were removed from the analysis. Moreover, only the half-monthly NDVI from January to September were used to derive the phenology indices (i.e. start and end of growing season length (GSL), see section 2.4).
The sum of monthly NDVI from January to December in a certain year was regard as the yearly NDVI (Gonsamo et al 2016(Gonsamo et al , 2017. Note that the illegitimate signals of the half-monthly NDVI data (i.e. NDVI<0.1) were filtered. Regional NDVI used in the trend analysis were averaged from the grids of the northern lands (>50°N). Before calculating the sensitivity of NDVI to MAT and the partial correlation between NDVI and seasonal temperature (see section 2.6), yearly NDVI data were resampled to raster of 0.5°×0.5°, to couple with the CRU temperature data.

Method of determining GSL
GSL was calculated as the difference between the start (SOS) and end (EOS) of growing season. The SOS and EOS were retrieved from the seasonal NDVI curve in each year based on the NDVI green-up thresholds, which were determined from the rate of seasonal changes in the multiyear mean NDVI (Piao et al 2006, Zhang et al 2013. More specifically, there were six steps in the determination of GSL. First, we calculated the seasonal curve of multiyear mean NDVI from 1982 to 2010 for each land grid cell and obtained the changing rate of NDVI (NDVI ratio ) as: where t is time throughout the year with an interval of 15 d. Then, after removing evident noise in the multiyear mean time-series curve of NDVI for each land grid cell, we performed a least-square regression analysis on the curves from January to September and from July to December for determining the NDVI thresholds of SOS and EOS, respectively, with an inverted parabola equation: where x is the Julian days and n is 6. The corresponding NDVI(t) with the maximum or minimum NDVI ratio was used as the NDVI threshold for determining SOS or EOS, respectively. Next, we performed a leastsquare regression analysis on the NDVI time-series curve in each year for each pixel after removing noise in the two different periods. After that, the SOS and EOS were identified from the fitted NDVI seasonal curves and their NDVI thresholds, by selecting the day when the fitted NDVI curve first reached the NDVI threshold. Finally, the GSL in each year for each land cell was calculated from the difference between EOS and SOS. The method in this study had been validated by ground-based phenological data in Tibetan Plateau  . The temporal anomalies were used for the linear-trend analyses. This analysis can provide both trend and its level of significance (i.e. the P value that quantifies the probability of whether the trend is statistically significant from zero) for each period. The moving-window method was used to detect whether the increasing trends of CO 2 indices are persistent. Comparing with the piecewise linear fitting method, it less depends on the results of single linear segment and the interval-length (Schleip et al 2008). This method has been used in detecting the changes of growing-season length and its response to climate change on various time scales (Rutishauser et al 2007, Schleip et al 2008, Jeong et al 2011, Fu et al 2015. Because the results based the moving-window analysis may be affected by the window-length (Fu et al 2015), we repeated the moving-window analyses with ten-, 15-and 20-year lengths.
The temporal trends of MAT (D MAT ), the sensitivity of NDVI to temperature (g NDVI ), and the sensitivity of CO 2 amplitude to temperature (g [ ] CO 2 amplitude ) in the periods of 1982-2010 and 1993-2007 were calculated in each grid cell. The g NDVI and g [ ] CO 2 amplitude were derived from the slope of linear regressions, representing the changes of NDVI and [CO 2 ] amplitude with per degree change of MAT. One-way ANOVA was used to estimate their differences between 1982-2010 and 1993-2007. All the analyses were applied in R (http://r-project.org/).
Partial correlation analyses were applied to exclude the impacts of the co-varying factors. For example, in calculating the impact of spring temperature on SOS, spring precipitation, spring solar radiation and last-year's autumn temperature were set as the controlling variables. Autumn precipitation, autumn solar radiation and spring temperature were set as the controlling variables to quantify the impact of autumn temperature on EOS. Similar method was used to detect the impact of spring and autumn temperature on annual NDVI and modeled GPP by replacing the seasonal precipitation and solar radiation with annual values. All the data were aggregated to the 0.5°× 0.5°resolution. The precipitation data were derived from the CRU TS4.0 dataset (Harris et al 2014) and the radiation data was from the Terrestrial Hydrology Research Group at Princeton University (Sheffield et al 2006).

Results and discussion
3.1. The temporal changes of atmospheric CO 2 seasonal cycle and vegetation greenness We first examined the trends of the measured annual CO 2 amplitude ([CO 2 ] amplitude ) at Point BRW, Alaska (71°N). The non-parametric Theil-Sen estimator showed that the increasing trends of the [CO 2 ] amplitude at the BRW (0.075 ppm yr −1 , P<0.05; figure 1(a)) were associated with the decreasing [CO 2 ] min (−0.058 ppm yr −1 , P<0.05) rather than the enhancing [CO 2 ] max (0.016 ppm yr −1 , P=0.17) from 1982 to 2010. The ten-year moving windows show that the increasing rates of the CO 2 amplitude was slower around 2000 (figure 1(a) and table S2). To avoid the biases from different time-intervals for trend estimation, we also detected the trends with 15-year (figure S3) and 20-year (figure S4) moving windows. The results also showed that the trends of [CO 2 ] amplitude from mid-1990 to mid-2000 (e.g. 0.03 ppm yr −1 , P=0.55, 1993-2007) (figure S4(a) and table S4) were significantly slower than those during other periods. A recent study which integrated the CO 2 records from multiple sites also has showed a stalled trend in the seasonality of atmospheric CO 2 during the same period (Yuan et al 2018).
A slowdown of vegetation greening since mid-1990s was also observed by our analysis on the dynamics of NDVI (figures 1(b) and S5(d)). This finding is consistent with the results from recent analyses on vegetation dynamics over the pan-Arctic tundra (Bhatt et al 2013, Jeong et al 2013). Both the MTE GPP (figure S5(e)) and the ground-based measurements of growing-season net ecosystem CO 2 exchange (figure S5(f)), the net ecosystem exchange data were derived from Belshe et al (2013) showed similar trends since 1990s. These lines of evidences together suggest that the increasing trend of the growing-season CO 2 uptake has weakened from 1990s to mid-2000s in northern ecosystems.
The synchronous changes of [CO 2 ] amplitude with [CO 2 ] min (figures S5(a)-(c)) and NDVI (figure 1) imply that the long-term positive trend of [CO 2 ] amplitude is, at least in part, driven by photosynthetic CO 2 uptake or vegetation growth (Forkel et al 2016, Wenzel et al 2016.  -1990s to mid-2000s (e.g. 1993-2007) than the whole period of 1982-2010, the sensitivities of [CO 2 ] amplitude and NDVI to MAT were lower during that period than 1982-2010. Given the fact that temperature in different seasons has non-uniform impacts on vegetation growth (Xia et al 2014), we further analyzed the changes of seasonal temperatures based on the CRU temperature datasets (see methods). As shown by figure 3(b), the fastest warming season was spring during 1985-1999 (+0.12°C year −1 ) but then changed to autumn during mid-1990s to mid-2000s (e.g. +0.11°C year −1 in 1993-2007, table S4). It indicates that a better understanding of the relationship between the seasonal temperature changes and vegetation growth is needed to explain the slowdown of [CO 2 ] amplitude from 1990s to mid-2000s.

The asymmetric responses of vegetation growth to spring and autumn warming
The variation of NVDI during 1982-2010 in northern ecosystems depends substantially on the GSL on both grid and regional scales (figures 1(b) and S6). Further partial correlation analysis showed that the SOS (partial r = −0.36) was more dependent on temperature change than the EOS (partial r = 0.018, figures 3(c) and S7). During 1982-2010, the SOS was advanced by 2.15 d°C −1 with spring warming, whereas warming in autumn only delayed EOS by 0.80 d°C −1 in northern ecosystems. As a result, the advancing rate of SOS over the moving 15 years has decreased but the extending rate of EOS was not significantly increased during 1982-2010 ( figure S8 and table S5). Meanwhile, the increasieng rate of NDVI until it stalled in mid-1990s is driven by warming-induced increase in spring and early summer NDVI along with the advancement of SOS (figure S8), for which spring warming has stalled after mid-1990s. It has contributed to the decline in the rate of [CO 2 ] amplitude increase since the mid-1990s. These results together suggest that the non-uniform warming between spring and autumn during mid-1990s and mid-2000s ( figure 3(b))   could be an important driving factor for the slowdown of expanding GSL, greening vegetation and the decreasing [CO 2 ] min . This can qualitatively explain the pause in the enhancement of [CO 2 ] amplitude ( figure 1(a)).

The response of vegetation productivity to spring and autumn warming in current TBMs
We examined whether TBMs that have been focusing on simulation of C dynamics in northern latitudes can adequately represent the differential impacts of spring and autumn warming on vegetation productivity. The ensemble output of GPP from five TBMs (CLM4.5, CoLM, ORCHIDEE, TEM6 and UVic; table S1) and a flux-based GPP dataset (MTE) were analyzed (figure 4). The dependence of modeled GPP variations on spring-temperature change (with the inter-model mean partial r as 0.57 under the significant level of P<0.05) is comparable with that of the MTE GPP (partial r=0.53, P<0.05) as well as that of NDVI (partial r=0.49, P<0.05; figure 4(b) and figure S9). However, the dependences of GPP variations on autumn-temperature change is more positive in the models (inter-model mean partial r=0.31, P < 0.05) than the MTE GPP (partial r=−0.04, P<0.05) and NDVI (partial r=−0.16, P<0.05; figures 4(b) and S10). This mismatch between modeling and dataoriented results indicates that the current TBMs overestimate the positive impact of rising MAT on ecosystem CO 2 uptake in the autumn.
As the land biogeochemical component in most Earth system models is similar to the TBMs in this study, it is still challenging to accurately simulate the seasonal cycle of atmospheric CO 2 . The findings of this study suggest that a better representation of the warming impacts on autumn phenology could partially improve the models' performance. However, the autumn phenology is diversely represented in different models. For example, leaf senescence in the ORCHIDEE model is simulated as the timing when monthly temperature falls below a given number, which varies with plant function type (Krinner et al 2005). In the TEM, growing season ended when the soil temperature is lower than the frozen point. However, leaf-senescence events are collectively affected by not only temperature but also day length (Ballantyne et al 2017), radiation (Bauerle et al 2012) and even spring phenology (Keenan and Richardson 2015, Liu et al 2016. In fact, the poor representation of autumn phenology by the models has been raised in some previous studies (Richardson et al 2010, 2012, Keenan and Richardson 2015. Thus, combining the different types of phenological data (e.g. Richardson et al 2018) with better phenology models could be helpful to improve the simulation of the seasonal cycle of atmospheric CO 2 at high latitudes in Earth system models.
3.4. The role of the non-uniform warming in regulating the seasonal atmospheric CO 2 cycle This study highlights that the asymmetric responses of vegetation growth to spring and autumn warming is an important driver for the decadal changes in the seasonality of atmospheric CO 2 . In spring, solar radiation is not limiting as temperature (Tanja et al 2003). Thus, spring warming extends the GSL by advancing the onset of plant photosynthesis (Piao et al 2008), leading to the increasing vegetation productivity and the decreasing [CO 2 ] min of the atmospheric CO 2 seasonal cycle. In autumn, solar radiation can Mean partial correlation coefficient (partial r) of NDVI and GPP to spring temperature and autumn temperature changes. The shaded areas represent the standard deviations of the partial r among models. Note that only the grid cells with significant partial correlation (P<0.05) were included in this analysis. obstruct the accumulation of abscisic acid (Gepstein and Thimann 1980) and substantially delay the timing of leaf senescence. Photoperiod is a more proximal factor than temperature in controlling senescence (Bauerle et al 2012). Thus, autumn warming has a neutral impact on vegetation productivity (figure 3(c)) over the northern lands.
Warming in autumn as well as in spring could potentially enhance the peak of atmospheric CO 2 seasonal cycle by stimulating the respiratory processes of plants and soil microorganisms (Piao et al 2008(Piao et al , 2017a. As shown by the FLUXCOM database (Jung et al 2017), the increasing trends of total ecosystem respiration during both growing and non-growing seasons were significantly larger in mid-1990s to mid-2000s (e.g. 1993-2007)

Conclusions
This study detected a slowdown of the increase in atmospheric CO 2 amplitude during mid-1990s to mid-2000s. This phenomenon was correlated with the pause of increasing NDVI and advancing SOS across the lands at northern high latitudes during the same period. The changes of vegetation greenness and growing-season length were temporally correlated with the stalled increase in spring temperature since mid-1990s. Warming in autumn was persistent during this period, suggesting that the non-growing season respiration could be more important in governing the future increase in seasonal CO 2 amplitude (Jeong et al 2018). These findings emphasize that the asymmetric responses of vegetation growth to spring and autumn warming is important in influencing the change of atmospheric CO 2 amplitude. This study also indicates that global carbon-cycle models need to better represent the phenological response to temperature change for accurately simulating the seasonal cycle of atmospheric CO 2 . Overall, this study confirms that the recent non-uniform climate warming among seasons has played an important role in regulating the temporal trends of vegetation growth and atmospheric CO 2 amplification over the northern lands.